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ABSTRACT 

1> 

The  acoustic  properties  of  marine  sediments  have  a  direct  effect  on  the 
propagation  of  sound  in  the  ocean.  In  the  frequency  range  of  interest  (50  -  500  Hz) 
the  sediment  can  be  modelled  as  a  fluid.  Assuming  horizontal  stratification  of  the 
ocean  bottom,  the  acoustic  parameters  of  interest  are  the  compressional  wave 
speed,  the  compressional  wave  attenuation  and  density  as  a  function  of  depth.  /  ^ >Jt 

P  r-r  J  S  -s 

SXn  inverse  method^  based  on  a  perturbation  technique^  presented  in  this^ 

r  ni  i » I  h-ij 

rjtiresra*  for  the  -  determination  ofy  these  parameters.  A  monochromatic  source 
experiment  is  proposed  because  of  the  desirability  of  such  an  experiment  for 
determining  the  acoustic  properties  of  an  anelastic  medium.  The  input  information 
is  the  plane  wave  reflection  coefficent  as  a  function  of  the  angle  of  incidence  at  a 
fixed  frequency.  A  nonlinear  integral  equation  relating  the  variations  of  these 
acoustic  properties  from  a  known  reference  value  to  the  plane  wave  reflection 

i)  f 

coefficient  is  derived,  -Tbis^is  then  linearised  using  the  Born  approximation.  The 
region  of  validity  of  the  Born  approximation  is  derived  and^  based  on  this^the 
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— 'optimum  angular  aperture  for  the  iuput  data  is  obtained. 

The  linearised  integral  equation  is  a  Fredholm  integral  equation  of  the  first 
kind.  An  acceptable  stable  solution  of  the  integral  equation  is  obtained  by  imposing 
a  priori  constraints  on  the  solution.  The  inversion  method  is  tested  using  synthetic 
data  and  inversions  are  carried  out  for  various  examples  of  the  attenuation 
coefficient  profile  and  the  sound  speed  profile.  The  results  obtained  with  noise  free 
data  show  good  agreement  between  the  true  profiles  and  the  reconstructed  profiles. 
The  resolution  obtainable  with  the  data  set  is  studied  using  the  resolving  power 
theory  of  Backus  and  Gilbert  and  the  inversion  method  is  shown  to  provide 
adequate  resolution.  The  effect  of  additive  noise  in  data  is  examined  and  inversions 
performed  with  noisy  data  yielded  stable  acceptable  results. 

Thesis  supervisor:  George  V.  Frisk,  Associate  Scientist,  Woods  Hole 
Oceanographic  Institution. 
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We  then  formulate  the  inverse  problem  for  the  determination  of  the 
attenuation  coefficient  profile  in  the  bottom.  This  is  done  in  Chapter  5.  Two 
approaches  to  linearisation  ,  namely  the  Born  and  Rytov  methods,  are  used  and  in 
each  case  an  integral  equation  relating  the  attenuation  coefficient  profile  to  the 
reflection  coefficient  is  derived. 

Chapter  6  is  a  predominantly  review  chapter,  where  we  study  the  problems  of 
instability  and  non-uniqueness  encountered  in  solving  inverse  problems  of  this  type. 
Various  methods  described  in  the  literature  for  solving  these  problems  are  then 
described  with  a  view  to  bringing  out  the  underlying  commonality  of  all  the 
approaches.  The  method  of  solving  nonlinear  problems  by  linearisation  is  then 
described. 

In  Chapter  7  we  continue  with  the  solution  of  the  inverse  problem  formulated 
earlier.  Adopting  a  regularisation  scheme  for  overcoming  problems  of  non¬ 
uniqueness  and  instability  several  examples  of  inversion  are  performed  using 
synthetic  data.  The  sensitivity  of  the  reconstruction  to  various  parameters  and 
other  related  issues  are  examined.  We  extend  the  formulation  of  the  inverse 
problem  to  determine  perturbations  in  other  acoustic  parameters  and  demonstrate 
this  by  simultaneously  solving  for  the  attenuation  profile  and  perturbations  in  the 
sound  speed  profile.  We  also  show  how  the  procedure  described  by  Backus  and 
Gilbert[20]  for  determining  the  resolvability  of  a  finite  data  set  can  be  applied  to 
situations  when  all  the  quantities  are  complex.  The  effect  of  noisy  data  on 
reconstruction  is  also  investigated. 

We  then  show  that  in  the  shallow  water  context,  the  problem  can  be 
reformulated  to  relate  the  perturbations  in  the  acoustic  parameters  to  the  changes 
in  the  eigenvalues.  This  is  done  in  Chapter  8.  Some  preliminary  results  of  inversion 
using  synthetic  and  real  data  are  presented. 
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perturbations  in  the  compressional  wave  speed  and  density.  We  do  this  not  only 
because  it  is  the  path  that  our  investigation  took  us  but  also  because  the  need  to 
determine  the  attenuation  coefficient  and  the  problems  in  its  determination  will 
then  be  brought  into  focus. 

When  the  ocean  is  a  shallow  water  wave  guide,  we  use  eigenvalues  of  the 
modes  trapped  in  the  water  column  as  input  information  for  the  inverse  method.  A 
linear  Fredholm  integral  equation  of  the  first  kind  relating  the  variations  in  the 
acoustic  parameters  from  a  known  reference  value  to  the  modal  eigenvalues  is 
obtained.  This  equation  can  then  be  solved  employing  any  one  of  the  methods 
available  for  solving  this  class  of  integral  equations. 

1.2  Overview 


s 


M 

$ 


■J 

i  ' 


$ 


In  Chapter  2  v/e  study  the  propagation  of  plane  waves  in  an  anelastic  medium 
to  obtain  a  correct  model  for  the  sediments.  This  is  important  as  it  can  influence 
the  type  of  experiment  performed  for  obtaining  the  acoustic  parameters. 

In  Chapter  3,  after  a  brief  review  of  the  methods  described  in  the  literature 
for  the  determination  of  the  attenuation  profile  in  marine  sediments,  we  propose  a 
method  where,  as  a  first  step  towards  inferring  the  acoustic  properties  of  the  ocean 
bottom,  we  obtain  the  plane  wave  reflection  coefficient  for  the  bottom  as  function 
of  angle  of  incidence[19].  A  sequence  of  steps  which  can  then  be  used  to  obtain  all 
the  acoustic  parameters  of  the  bottom  is  indicated. 
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In  Chapter  4  we  describe  methods  for  solving  the  forward  problem  i.e  finding 
the  acoustic  field  in  the  bottom  given  its  acoustic  parameters.  A  numerically  stable 
propagator  matrix  algorithm  is  presented. 
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distinguish  it  from  other  Boro  approximation  methods  are 

1.  The  input  information  is  the  plane  wave  reflection  coefficient  as  a 
function  of  the  angle  of  incidence  at  a  fixed  frequency. 

2.  The  reference  wave  speed,  density  and  the  wave  attenuation  can  vary 
arbitrarily  with  depth. 

3  The  medium  can  be  lossy. 

4.  The  method  yields  all  three  acoustic  parameters. 

5.  The  Born  approximation  is  not  uniformly  valid  over  all  the  angles  and 
therefore  band  limited  information  is  used. 

The  linear  integral  equation  obtained  by  this  method  is  a  Fredholm  integral 
equation  of  the  first  kind.  These  equations  suffer  from  problems  of  non-uniqueness 
and  instability  of  the  solution.  We  obtain  a  stable  acceptable  solution  by  using  the 
regularisation  method  due  to  Phillips[l7]  and  Twomey[18]. 

Though  the  problem  as  described  above  involves  the  determination  of  all 
three  acoustic  parameters,  the  starting  point  of  this  thesis  was  the  development  of 
a  method  for  obtaining  the  compressional  wave  attenuation  for  the  ocean  bottom. 
In  studying  a  possible  method,  we  noted  that  it  required  an  exact  knowledge  of  the 
sound  speed  and  density  profiles  for  the  bottom.  We  found  that  even  small  errors  in 
these  two  parameters  adversely  affected  the  determination  of  the  attenuation 
coefficient.  This  led  to  the  present  approach  where  instead  of  trying  to  correct  for 
the  errors  in  the  sound  speed  and  density,  these  are  treated  as  unknowns  in  a 
general  formulation  that  yields  corrections  to  the  sound  speed  and  density  profile* 
together  with  the  attenuation  coefficient  profile. 

We  will,  therefore,  develop  the  method  for  the  determination  of  the 
attenuation  coefficient  profile  for  the  bottom  and  then  extend  it  to  deal  with 
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through  a  linear  integral  equation  as  before. 

Subsequent  to  the  publication  of  the  first  paper  by  Cohen  and  Bleistein,  a 
number  of  papers  have  appeared  in  the  field  of  seismology  which  are  based  on  a 
Born  type  approximation(5-9j.  Features  common  to  all  these  are  as  follows. 

1.  The  response  to  an  impulsive  source  at  the  surface  of  the  medium 
being  probed  is  measured  at  the  same  location  as  the  source  or  at  a 
fixed  off-set  from  it.  This  is  equivalent  to  measuring  the  plane  wave 
reflection  coefficient  as  a  function  of  frequency  at  normal  incidence  or 
at  a  fixed  angle  of  incidence. 

2.  The  reference  wave  speed  and  density  are  generally  assumed  to  be 
constants.  Raz[10],Clayton  and  Stolt[ll],  and  Bleistein  and  Gray[12], 
however,  also  consider  the  case  when  the  reference  compressional  wave 
speed  is  depth  dependent. 

3.  The  medium  is  lossless. 

4.  The  variation  of  true  wave  speed  from  the  reference  speed  is  small  and 
the  Born  approximation  is  valid  at  all  frequencies  of  the  broad  band 
source. 

The  experiment  that  we  propose  uses  a  CW  source  and  the  input  data  in  this 
case  is  the  plane  wave  reflection  coefficient  at  various  angles  of  incidence.  A 
schematic  of  the  experiment  is  shown  in  Figure  1-3.  Exact  methods  for  obtaining 
the  acoustic  parameters  of  the  medium  from  the  angular  dependence  of  the 
reflection  coefficient  have  been  proposed  by  Hooshayar  and  Razavy[l3]  and 
Stickler[14|.  Schaubert  and  Mittra[15j  and  Roger[16]  present  inverse  methods  based 
on  perturbation  techniques  for  determining  the  permittivity  of  a  lossless  dielectric 
medium  using  a  monochromatic  source  to  probe  the  medium.  However,  a 
perturbation  method  using  the  angular  dependence  of  the  plane  wave  reflection 
coefficient  as  input  information  has  not  been  applied  to  the  determination  of  all 
three  acoustic  parameters  of  the  ocean  bottom.  The  features  of  this  method  that 


have  a  direct  effect  on  the  propagation  of  sound  in  the  ocean. 

The  forward  scattering  problem  deals  with  the  determination  of  the  wave 
field  scattered  from  an  object  with  known  properties  for  a  given  incident  wave 
field.  In  contrast  to  this,  in  the  inverse  scattering  problem  the  scattered  field  for  a 
given  incident  field  is  measured  at  the  boundary  of  the  object  or  at  points  external 
to  the  object  and  from  this  the  properties  of  the  object  are  inferred.  The 
determination  of  the  acoustic  properties  of  the  ocean  bottom  from  the  plane  wave 
reflection  coefficient  falls  into  this  category. 

Two  approaches  are  found  in  the  literature  for  determining  the  acoustic 
properties  of  a  layered  medium  from  its  reflection  response.  These  are  exact 
methods  and  approximate  methods.  A  review  of  exact  inversion  methods  using 
reflection  coefficient  data  for  layered  media  is  given  by  Newton[2].  Our  approach  is, 
however,  based  on  an  approximate  method  developed  by  Cohen  and  Bleistein[3j.  In 
the  setting  of  geophysical  inverse  problems,  their  method  was  applied  to  the 
determination  of  the  compressional  wave  speed  in  the  interior  of  the  earth.  A 
schematic  diagram  for  the  experiment  is  shown  in  Figure  1-2.  The  source  is  an 
impulsive  plane  wave  source.  The  backscatterd  field  at  the  location  of  the  source  is 
measured.  Using  perturbation  theory,  the  back  scattered  field  is  related  to  the 
variation  of  the  compressional  wave  speed  from  a  known  reference  value  and  a 
linear  integral  equation  obtained.  The  linearising  approximation  is  called  the  Born 
approximation.  When  the  reference  wave  speed  is  a  constant  a  closed  form  solution 
to  the  integral  equation  is  obtained.  In  the  two  dimensional  case[4],  the  experiment 
consists  of  using  a  line  of  sources  which  is  equivalent  to  repeating  the  experiment 
for  the  one  dimensional  problem  at  a  number  of  locations  along  a  line.  In  this  case 
the  variation  of  the  compressional  wave  speed  from  its  reference  is  in  two  spatial 
coordinates.  This  is  related  to  the  backscattered  field  measured  along  the  line 
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Chapter  1 
Introduction 


The  prediction  of  sound  propagation  in  the  ocean  requires  a  knowledge  of  the 
acoustic  properties  of  marine  sediments.  In  the  range  of  frequencies  of  interest  (50 
-  500  Hz)  marine  sediments  can  be  modelled  as  a  fluid[lj.  Therefore,  the  acoustic 
properties  of  marine  sediments  of  interest  are  the  compressional  wave  speed,  the 
compressional  wave  attenuation  and  the  density.  Further  when  the  assumption  of 
horizontal  stratification  is  made,  these  parameters  are  a  function  of  depth  only.  In 
this  thesis,  we  present  an  inverse  method  based  on  a  perturbation  technique  for 
obtaining  these  parameters.  The  input  information  is  the  plane  wave  reflection 
coefficient  for  the  bottom  as  a  function  of  angle  of  incidence  at  a  fixed  frequency. 


1.1  Background 

Figure  1-1  is  a  typical  ray  diagram  for  low  frequency  sound  (50  -  500  Hz) 
interacting  with  the  bottom.  The  rays  emanate  from  the  source  and  are  partially 
reflected  at  the  water/ocean  bottom  interface  and  the  energy  entering  the  bottom, 
therefore,  depends  on  the  impedance  contrast  at  this  interface.  The  rays  that 
penetrate  into  the  bottom  are, however,  refracted  back  into  the  water  column 
because  of  the  increasing  sound  speed.  As  the  rays  travel  in  this  region,  some  of  the 
energy  associated  with  them  is  converted  into  heat.  The  attenuation  coefficient  is  a 
measure  of  the  amount  of  energy  dissipated  in  this  manner.  Thus  the  compressional 
wave  speed,  the  compressional  wave  attenuation  and  the  density  of  the  sediment 


•Tjufc  TW  T  .*■**.  ‘ 


List  of  Tables 


Table  5-1:  Real  and  imaginary  parts  of  the  wavenumber  at  220  Hz  for  9-4 
marine  sediments 

Table  7-1:  Parameters  used  for  ocean  bottom  model  164 

Table  7-II:  Resolution  lengths  for  the  different  examples  of  189 
reconstruction 

Table  7-III:  Resolution  length  and  aperture  size  191 

Table  7-IV:  Resolution  length  and  number  of  data  points  192 

Table  7-V:  Exact  and  erroneous  sound  speed  profile  200 

Table  7-V1:  True  and  background  model  parameters  -  case(i)  204 

Table  7-VII:  True  and  background  model  parameters  -  ease(ii)  205 


Figure  7-35: 
Figure  7-36: 
Figure  7-37: 
Figure  7-38: 
Figure  7-39: 
Figure  7-40: 
Figure  7-41: 

data-40dB 
Figure  7-42: 
-40dB 

Figure  7-43: 
40dB 

Figure  7-44: 
50dB 

Figure  7-45: 
50dB 

Figure  7-46: 
50dB 


Effect  of  error  in  sound  speed  on  reconstruction 
Effect  of  error  in  density  on  reconstruction 
Simultaneous  reconstruction  •  sound  speed  profile;  case  (i) 
Simultaneous  reconstruction  -  Attenuation  profile;  case(i) 
Simultaneous  reconstruction  -  Sound  speed  profile;  case(ii) 
Simultaneous  reconstruction  -  Attenuation  profile;  case(ii) 
Reconstruction  of  profile  in  figure  7-1 4(a)  with  noisy 

Reconstruction  of  profile  in  figure  7-14(b)  with  noisy  data 

Reconstruction  of  profile  in  figure  7-14(c)  with  noisy  data- 

Reconstruction  of  profile  in  figure  7*14(a)  with  noisy  data 

Reconstruction  of  profile  in  figure  7-14(b)  with  noisy  data 

Reconstruction  of  profile  in  figure  7- 14(c)  with  noisy  data 


Figure  8-1: 
Figure  8-2: 
Figure  8-3: 
Figure  8-4: 
Figure  8-5: 
Figure  8-6: 
Figure  8-7: 
Figure  8-8: 

case(i) 
Figure  8-9: 

case(ii) 
Figure  8-10: 


Experimental  configuration  for  shallow  water  experiment 
Hankel  transform  of  measured  pressure  field 
Ocean  model 

Exact  and  guess  model  -  casefi) 

Exact  and  guess  model  -  case(ii) 

Reconstructed  profile  -  casefi) 

Reconstructed  profile  -  case(ii) 

Reconstructed  profile  using  data  from  all  three  frequencies- 
Reconstructed  profile  using  data  from  all  three  frequencies- 
Inversion  of  real  data 


198 

199 
206 

207 

208 
209 
211 

212 

213 

214 

215 

216 

222 

223 

225 

230 

231 

232 

233 

234 

235 


237 


Chapter  2 

Acoustic  Modelling  of  Marine  Sediments 


Acoustic  transmission  in  the  ocean  at  low  frequency  is  dependent  on  the 
geoacoustic  properties  of  the  seabed,  namely,  the  speed  and  attenuation  of  the 
compressional  and  shear  waves  and  density.  A  complete  specification  of  the  seabed 
will  contain  the  wave  speeds  and  attenuations  as  a  function  of  the  three  spatial  co¬ 
ordinates  and  frequency,  and  density  as  a  function  of  the  spatial  co-ordinates. 

In  this  chapter  we  develop  a  model  for  the  ocean  bottom  that  will  be  used  in 
the  analysis  to  be  presented  in  the  following  chapters  of  this  work.  In  this  context, 
we  study  the  frequency  dependence  of  attenuation  and  the  dispersion  of  the  wave 
speeds  due  to  attenuation.  These  are  important  considerations  in  the  design  of  a 
suitable  experiment  for  the  measurement  of  the  attenuation  coefficient.  We  also 
summarise  the  other  geoacoustic  characteristics  associated  with  wave  propagation 
in  marine  sediments. 

2.1  Dispersion  and  attenuation  of  compressional  waves 

We  start  by  studying  the  propagation  of  plane  compressional  waves  in  an 
attenuating  medium.  In  deriving  the  classical  loss-less  wave  equation  it  is  assumed 
that  the  material  behaves  in  an  elastic  manner.  The  result  of  such  an  assumption  is 
that  sound  waves  propagate  in  the  medium  without  change  of  shape  and  suffer  no 
attenuation.  However,  in  real  sediments  the  waves  decay  as  they  propagate.  Their 
energy  is  gradually  converted  to  heat.  This  process  is  the  result  of  a  variety  of 
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mechanisms  which  operate  for  most  part  at  the  microscopic  level.  At  the 
macroscopic  level  the  energy  loss  is  given  such  terms  as  internal  friction, 
attenuation  and  anelasticity. 

The  most  direct  method  of  defining  ’internal  friction’  is  as  the  ratio  —(AEJ)jE 
where  AE  is  the  energy  lost  in  each  cycle  because  the  material  is  not  perfectly 
elastic  and  E  is  the  peak  energy  stored  in  the  volume  of  material  taken  through  the 
stress  cycle  at  a  frequency  u.  The  specific  attenuation  coefficient  Q  is  defined  as 

2ff  AE 

—  = -  (2.1) 

Q  E 

This  factor  Q  can  be  measured  without  any  assumption  being  made  on  the 
attenuating  mechanisms  and  is  therefore  used  to  describe  attenuating  behaviour  of 
the  medium  in  a  phenomenological  sense. 

For  a  plan?  propagating  wave,  the  effect  of  ’internal  friction’  will  be  the 
gradual  decay  of  the  amplitude  of  the  wave  as  it  propagates.  Since  the  energy  is 
propotional  to  the  square  of  the  amplitude, 


Q(uj)  2irE  ir A 

In  the  above  expression  AA  is  the  change  in  amplitude  over  one  cycle  or 
AA=(dA/dx)\  where  X  is  the  wave  length. 


dA  2irC(u) 

AA  = - 

dx  u 


A(x)  =  \exp( - ) 

^  2  C(u)Q(u) 

Therefore,  for  a  propagating  wave, 


rtCQUCMCT.  kM* 


Figure  2-1:  Variation  of  the  attenuation  coefficient  with 
frequency-reproduced  from  Hamilton[4] 
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A(x,t)  =  A^exp[i{——+i  --  ---  -wO] 
C[v)  Q(<jj)2C(<jj) 

=Aq  exp  [  i  (*(w)-wf] 


Here  k(u)  is  a  complex  quantity,  its  imaginary  part  representing  attenuation. 
Another  measure  of  attenuation  found  in  the  literature  is  the  attenuation 
coefficient  which  is  given  by  the  following  expression. 


a(cj)  = 


2C((j)Q(u) 


Attenuation  in  a  medium  can  therefore  be  accounted  for  by  making  the  wave 
number  complex. 

The  attenuation  coefficient  is  a  function  of  frequency  by  definition.  Aki  and 
Richards[l]  show  that  the  behaviour  of  waves  propagating  in  an  attenuating 
medium  cannot  be  explained  unless  the  assumption  of  dispersion  of  these  waves  is 
made.  Therefore  both  Q  and  C  are  functions  of  frequency.  Figure  2-1  is  a  plot  of 
the  attenuation  coefficient  for  compressional  waves  with  respect  to  frequency  based 
on  laboratory /field  experiments.  Most  of  the  information  is  in  the  high  frequency 
region.  The  experimental  results  show  a  rough  linear  dependence  of  the  attenuation 
coefficient  on  frequency  in  the  high  frequency  range. 

Efforts  to  explain  this  experimentally  observed  behaviour  have  been  made  by 
modifying  the  equations  of  th  classical  elastic  theory.  Different  approaches  have 
been  proposed.  For  example  Hamilton[4]  proposes  a  visco-elastic  model  for  water 
saturated  sediments  while  another  approach,  called  the  Biot  model,  which  treats 
the  sediment  as  a  porous  elastic  solid  saturated  by  a  viscous  fluid  has  been  pursued 
by  Stoll{13].  The  two  models  will  be  described  briefly  in  the  following  paragraphs. 


2.1.1  Visco-elastic  model 


In  the  classical  elastic  theory,  the  stress  and  the  strain  are  related  by  a  linear 
relationship  such  as 

o  ^Be  (2.7) 

where  a  is  the  stress,  c  the  strain  and  B  is  the  constant  of  propotionality  called  the 
elastic  constant.  In  the  visco-elastic  model  the  anelasticity  is  incorporated  by 
modifying  the  elastic  equation  to  express  the  stress  as  a  function  of  strain  and 
strain  rate  as  indicated  in  equation  (2.8),  where  dt/dt  is  the  strain  rate. 

dt 

a  =  Bt+G  —  (2.8) 

dt 

We  now  have  two  elastic  constants  B  and  G.  A  mechanical  system  representing  this 
stress  strain  relationship  is  shown  in  Figure  2-2  The  spring  represents  a  perfectly 
elastic  body  for  which  the  stress-strain  relationship  is  governed  by  the  law 

a  =Bt  (2.9) 

The  dashpot  represents  a  perfectly  viscous  body  for  which  the  stress  strain 
relationship  is  given  below. 

dt 

a  =  —  (2.10) 

dt 

For  a  system  in  which  these  two  elements  are  in  parallel  the  relationship  in 
equation  (2.8)  is  obtained.  The  analogies  of  various  other  visco-elastic  models  to 
mechanical  systems  are  given  in  Reference  5. 

The  model  in  Figure  2-2  is  called  the  Kelvin-Voigt  model  and  Hamilton[2,4] 
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has  proposed  this  model  to  describe  the  anelastic  behaviour  of  marine  sediments. 
Kolsky[3]  has  shown  that  the  elastic  equation  (2.8)  when  applied  to  solids  resulted 
in  the  Lame’s  constants  Xand  /i  being  replaced  by  X+X'  and  p+p'  where  X'  and  p1 
govern  the  energy  damping  and  X  and  p  together  with  the  density  govern  the 
velocity.  Based  on  this  analysis,  the  following  equations  for  specific  attenuation 
coefficient  and  compressional  wave  speed  were  obtained  by  Hamilton(2,4). 

1  _ 

Qif)  ~ 

o  (!-f2) 

(X+2/i)  =  pc  2  — — 
v  {i+ry 

where  r  =  a(J)C(f)/2nf.  When  damping  is  small  r  is  neglected  and  we  recover  the 
equation  (2.6).  Again  when  r  is  negligible  there  is  no  dispersion  in  the  medium.  For 
this  visco-elastic  model  the  attenuation  coefficient  is  proportional  to  the  square  of 
frequency  [3). 

The  Kelvin-Voigt  model  does  not  exhibit  the  behaviour  of  real  materials  when 
a  load  is  applied  to  it.  The  response  of  the  Kelvin-Voigt  model  to  a  suddenly 
applied  stress  is  shown  in  Figure  2-3.  The  strain  is  zero  at  t=0  and  reaches  its  final 
value  asymptotically  with  time.  However,  real  materials  show  instantaneous  strain. 
A  model  which  overcomes  this  shortcoming  is  the  "standard  linear  model”.  In  this 
model  the  stress  and  stress  rate  are  expressed  as  a  function  of  strain  and  strain  rate 
and  is  given  by  equation  (2.11) 

da  dt 

(211) 

where  MR  is  called  the  deformation  modulus,  r(  is  the  stress  relaxation  time  under 
constant  strain  and  r  is  the  strain  relaxation  time  under  constant  stress.  We  now 
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have  three  elastic  constants.  A  mechanical  system  that  models  this  behaviour  and 
its  response  to  a  suddenly  applied  load  are  in  Figures  2-4  and  2-5.lt  is  seen  from 
Figure  2-5  that  there  is  an  instantaneous  response  of  this  system  to  suddenly 
applied  load  as  is  the  case  for  real  materials.  The  dependece  of  the  specific 
attenuation  coefficient  for  this  model  is  shown  in  Figure  2-6.  The  Kelvin-  Voigt 
model  anu  the  standard  linear  model,  therefore,  cannot  explain  the  linear 
relation&nip  found  experimentally. 

However, since  energy  damping  is  due  to  a  large  number  of  mechanisms,  Liu 
etal[6]  proposed  that  by  considering  a  large  number  of  relatively  closely  spaced 
relaxation  mechanisms  representing  the  various  damping  processes,  with  each 
mechanism  behaving  like  a  linear  solid,  the  attenuation  coefficient  can  be  shown  to 
be  linearly  related  to  frequency  for  earth  like  materials. 

Using  the  theory  of  superposition  of  relaxation  mechanisms  Liu  etal[6]  derived 
the  following  expression  for  the  dispersion  of  waves. 


q,K> 


1 

1+ — 

*Q 


(2.12) 


A  similar  equation  for  dispersion  has  been  derived  by  other  investigators(8,9]  using 


entirely  different  approaches. 


The  theory  of  superposition  of  relaxation  mechanisms,  therefore,  gives  the 
experimentally  observed  linear  dependence  of  attenuation  with  frequency  and  also 
provides  a  dispersion  relationship  which  has  been  obtained  from  two  completely 
different  approaches.  This  theory  can,  therefore,  be  considered  as  providing  the 
physical  basis  for  understanding  the  propagation  of  waves  in  an  anelastic  medium 
like  the  earth.  Treating  consolidated  marine  sediments  as  earth-like  materials  and 
applying  the  above  theory,  the  attenuation  coefficient  can  be  predicted  to  have  a 


B 
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linear  relationship  with  respect  to  frequency  over  the  entire  frequency  range  and 
the  values  of  the  attenuation  coefficient  obtained  at  high  frequency  can  then  be 
extrapolated  into  the  low  frequency  region.  This  procedure  has  been  recommended 
by  Hamilton[4|.  Experimental  results  obtained  by  Frisk[18),  Mitchell  and  Focke[19] 
and  Stoll[20j  for  attenuation  in  marine  sediments  at  low  frequencies,  however,  yield 
values  much  less  than  what  is  predicted  by  extrapolating  the  linear  law.  This  casts 
doubt  on  the  applicability  of  the  superposition  theory  to  marine  sediments. 


2.1.2  ’  Biot  ’  model 


Biot[10,ll,12]  studied  the  propagation  of  elastic  waves  in  a  system  composed 
of  a  a  porous  elastic  solid  saturated  by  a  viscous  fluid.  The  fluid  is  assumed  to  be 
compressible  and  may  move  relative  to  the  frame.  In  this  model,  the  losses  are 
grouped  into  two  fundementally  different  categories,  one  which  accounts  for  the 
anelasticity  of  the  skeletal  frame  and  the  other  for  the  viscosity  of  the  fluid  moving 
relative  to  the  frame. 

The  theory  proposed  by  Biot  was  later  applied  by  Stoll[13]  and  Stoll  and 
Bryan[14]  for  propagation  in  ocean  sediments.  Using  this  theory,  the  variation  of 
compressional  wave  speed  with  frequency  and  the  variation  of  attenuation  with 
frequency  can  be  obtained  for  known  sediment  properties. 

The  way  in  which  attenuation  varies  with  frequency  depends  on  the 
dominance  of  one  or  the  other  of  the  two  modes  of  energy  dissipation  that  are  built 
into  the  model,  namely,  the  viscous  losses  in  the  fluid  as  it  moves  relative  to  the 
frame  or  the  frictional  losses  in  the  skeletal  frame.  Figure  2-7  shows  two  kinds  of 
response  predicted  by  the  model.  For  very  fine  materials  of  low  permeability,  such 
as  silty  clay,  the  losses  in  the  skeletal  frame  dominate  and  the  authors  propose  a 


linear  relationship  while  for  high  permeability  material  like  coarse  sand  etc.  a 
dependence  /“with  n  =  2  at  the  low  frequency  end  and  n  =  1/2  at  the  high 
frequncy  has  been  proposed.  McDaniel  and  Beebe[15j  applied  the  Biot  theory  and 
computed  the  acoustic  characteristics  of  different  types  of  sediments.  The 
sediments  were  classified  on  the  basis  of  their  mean  grain  size.  To  derive 
expressions  for  the  acoustic  characteristics  they  used  empirical  relations  between 
permeability  and  grain  size.  The  variation  of  attenuation  with  frequency  was  found 
to  obey  the  law  a(J)  —  Af1  where  n  lies  between  1  and  1.8.  The  results  obtained 
by  them  are  shown  in  Figure  2-8.  Hovem[l6],  however,  suggests  that  the  viscous 
losses  which  is  the  cause  of  the  nonlinear  dependence  can  be  made  to  vary  linearly 
with  frequency  if  a  distribution  of  grain  sizes  is  adopted. 

Figure  2-11  is  a  collection  of  laboratory  and  field  data  for  attenuation  plotted 
against  frequency  with  the  prediction  of  the  Biot  model  as  computed  by  Stoll[23] 
superimposed.  We  note  that  in  the  low  frequency  region  there  is  considerable 
scatter  in  the  experimentally  determined  values  and  applicability  of  the  Biot  model 
to  marine  sediments  cannot  be  confirmed  on  the  basis  of  the  evidence  presented  in 
this  figure. 

The  variation  of  compressional  wave  with  frequency  as  obtained  using  the 
Biot  model  is  shown  in  Figure  2-9  for  sediments  of  different  mean  grain  size.  For 
coarser  materials  (low  <j>)  the  dispersion  is  appreciable,  which  is  similar  to  the  result 
obtained  by  Liu  as  such  materials  also  have  low  Q. 

2.1.3  Compressional  wave  attenuation  vs  Depth 

The  ray  diagram  in  Figure  1-1  shows  the  rays  interacting  with  the  bottom. 
The  loss  of  energy  suffered  due  to  interaction  with  the  bottom  will,  therefore, 
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Marine  sediments  have  enough  rigidity  to  transmit  shear  waves.  When 
compressional  waves  are  incident  on  the  ocean  bottom  interface,  shear  conversion 
takes  place.  If  the  ocean  bottom  is  modelled  as  a  set  of  layers  with  velocity 
gradients  within  each  layer,  then  such  shear  conversion  will  occur  at  the  layer 
interfaces  and  if  the  velocity  gradients  are  large  it  will  occur  continuously  within 
the  layers.  Fryer[17]  studied  the  effect  of  shear  in  marine  sediments  and  concludes 
that  the  effect  of  shear  conversion  within  the  the  sediment  is  small  at  frequncies 
above  20  Hz.  However,  if  strong  discontinuities  in  shear  speeds  exist  or  if  any 
energy  reaches  the  basement  where  the  shear  speed  becomes  comparable  to  the 
compressional  wave  speed,  then  shear  conversion  cannot  be  ignored.  For  thick 
sediments,  we  can  assume  that  there  are  no  strong  discontinuities  in  the  shear 
speed  within  the  sediment  and  very  little  energy  reaches  the  basement.  Under  such 
assumptions,  the  effect  of  shear  conversion  can  be  ignored  and  the  sediment  treated 
as  a  fluid.  The  ocean  bottom  model  now  contains  only  three  parameters,  the 
compressional  wave  speed,  compressional  wave  attenuation  and  density. 

Some  investigators  have  made  a  further  simplifying  assumption  that  the 
dispersion  of  compressional  waves  is  negligible  and  therefore  the  compressional 
wave  speed  is  a  function  of  depth  alone.  For  reasons  explained  earlier  we  will  not 
make  such  an  assumption.  Further  we  do  not  assume  linear  dependence  of 
attenuation  with  frequency  since  the  applicability  of  this  behaviour  at  low 
frequencies  to  all  types  of  sediments  is  in  doubt. 

The  model  for  the  ocean  bottom  we  adopt  is  shown  in  figure2-6.  The  figure 
represents  a  thick  sediment  layer  in  which  the  compressional  wave  speed  ,  density 
and  attenuation  can  vary  arbitrarily  with  depth,  overlying  a  sub-bottom  of 
constant  density,  compressional  wave  speed  and  attenuation.  Since  we  assume  that 
very  little  energy  reaches  the  basement,  the  sub-bottom  is  modelled  as  an  infinite 
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Figure  2-17:  Ocean  bottom  model 
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Figure  2-15:  Variation  of  density  with  depth 
reproduced  from  Hamilton  [4] 


permeability  material.  In  the  context  of  the  evidence  available  for  compressional 
waves  the  validity  of  this  assumption  in  the  low  frequency  region  is  questionable. 
Since  the  attenuation  coefficient  for  shear  waves  is  at  least  an  order  of  magnitude 
larger  than  that  of  compressional  waves,  dispersion  in  this  case  will  be  substantial. 

2.2.4  Shear  wave  attenuation  as  a  function  of  depth 

Very  little  information  is  available  on  the  variation  of  shear  wave  attenuation 
with  depth.  For  modelling  purposes  Hamilton  recommends  that  the  shear  wave 
attenuation  be  varied  with  depth  in  proportion  to  the  variation  of  compressional 
wave  attenuation  with  depth. 

2.2.5  Density  as  a  function  of  depth 

Density  of  the  sediments  can  be  modelled  as  increasing  monotonically  with 
depth  due  to  overburden  pressure.  Variation  of  density  with  depth  is  in  Figure 
2-15 

2.3  Ocean  bottom  model 

In  a  complete  model  for  the  ocean  bottom  the  geoacoustic  parameters  are 
functions  of  the  spatial  co-ordinates.  A  simplifying  assumption  commonly  made  for 
the  ocean  bottom  model  is  that  of  horizontal  stratification  i.e.,  the  geoacoustic 
parameters  depend  spatially  only  on  depth.  Figure  2-16  is  a  section  of  the  ocean 
bottom  obtained  by  a  3.5  kHz  seismic  profiler  in  the  Icelandic  Basin[22).  The  well 
defined  layering  of  the  bottom  can  be  seen  in  this  figure  and  such  layering  is 
observed  in  many  regions  of  the  ocean.  Thus  the  assumption  of  horizontal 
straticication  is  frequently  satisfied. 


2.2  Summary  of  other  essential  characteristics  of  the  geoacoustic 
parameters. 

2.2.1  Compressional  wave  speed  as  a  function  of  depth 

In-situ  measurements  of  velocity  at  the  sediment  water  interface  show  a  small 
discontinuity  in  the  compressional  wave  speed.  The  velocity  ratio  defined  as  the 
ratio  of  sound  speed  in  the  sediment  to  the  sound  speed  in  water  range  from  0.97  to 
1.03.  From  the  sediment  water  interface  the  sound  speed  increases  monotonically. 
Using  experimental  data  on  compressional  wave  speed  at  various  depths 
Hamilton[4]  gives  a  regression  equation  for  modelling  the  compressional  wave  speed 
in  the  sediment.  The  variation  of  compressional  wave  speed  with  depth  is  plotted  in 
Figures  2-12and  2-13  for  different  types  of  sediments. 

2.2.2  Shear  wave  speed  as  a  function  of  depth 

All  marine  sediments  possess  enough  rigidity  to  transmit  shear  waves.  Shear 
waves  are  important  in  underwater  sound  propagation  because  compressional  waves 
can  be  partially  converted  to  shear  waves  at  reflection  boundaries  and  by 
compressional  and  shear  velocity  gradients.  Variation  of  shear  wave  velocity  with 
depth  is  given  in  Figure  2-14.  For  modelling  purposes  Hamilton[4)  gives  a  regression 
equation  that  relates  the  compressional  wave  speed  to  shear  wave  speed. 

2.2.3  Shear  wave  attenuation  and  dispersion 

Experimental  results  for  shear  wave  attenuation  are  far  fewer  than  that  of 
compressional  waves.  For  modelling  purposes  Hamilton[4]  suggests  the  use  of  the 
linear  relationship  between  attenuation  and  frequency  in  the  case  of  high 
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depend  on  the  sediment  layers  sampled  by  these  rays.  Since  each  ray  has  different 
ray  paths  the  energy  ioss  associated  with  it  will  depend  on  the  variation  of 
attenuation  with  depth.  Very  little  experimental  data  is  available  on  this  important 
characteristic.  The  available  data  is  shown  in  Figure  2-10  which  is  a  plot  of  fc(z). 
The  function  k(z)  is  derived  from  the  attenuation  coefficient  using  the  law 
a(f,z)=k(z)f.  In  the  case  of  silty-clay  sediments  a  peak  in  the  value  of  attenuation 
is  observed  around  300  m.  Recent  field  experiments  carried  out  by  Jacobson  etal[21] 
and  Mitchell  and  Focke[19j  also  indicate  the  existence  of  such  peaks. 


2.1.4  Conclusions 

The  modelling  of  the  sediment  as  a  visco-elastic  model  or  as  a  porous  solid 
leads  to  conflicting  results  especially  in  the  low  frequency  region.  Experimental 
results  in  the  low  frequency  region  is  limited  and  the  question  as  to  which  of  these 
models  is  applicable  for  marine  sediments  has  not  yet  been  resolved.  However, 
based  on  the  theory  presented  so  far,  the  following  conclusions  can  be  drawn. 


1.  In  the  high  frequency  region  (above  1  KHz),  attenuation  in  low 
permeability  materials  can  be  modelled  as  having  a  linear  relationship 
with  frequency. 

2.  The  behaviour  of  marine  sediments  in  the  low  frequency  region  is  not 
well  understood  and  under  such  circumstances  it  will  not  be  correct  to 
extrapolate  the  results  obtained  in  the  high  frequency  region  to  the  low 
frequencies. 

3.  Dispersion  of  waves  for  coarse  sediments  can  be  substantial  especially 
over  a  wide  range  of  frequencies. 
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Chapter  3 

Determination  of  the  Attenuation  Profile 


Attenuation  as  described  in  Chapter  2  refers  to  the  decrease  in  the  amplitude 
of  propagating  waves  due  to  a  variety  of  energy  damping  mechanisms.  This 
attenuation,  also  called  intrinsic  attenuation,  depends  on  the  type  of  sediment  (i.e 
the  compositon  of  the  sediment,  the  particle  size  and  other  parameters).  Laboratory 
measurements  of  intrinsic  attenuation  can  be  readily  performed  but  are  restricted 
to  high  frequencies  because  of  the  long  wave  lengths  at  low  frequencies  and  the 
difficulty  in  eliminating  boundary  effects.  Laboratory  measurements  have, 
therefore,  been  limited  to  frequencies  above  1  KHz  [1-5]. 

In  view  of  this,  field  experiments  which  utilize  the  in-situ  deposits  of 
sediments  as  large  scale  specimens  have  been  performed.  In  such  field  experiments 
it  is  necessary  to  take  into  account  all  factors  like  spreading,  reflections,  etc.  which 
affect  the  amplitude  of  the  propagating  waves.  The  determination  of  the 
attenuation  coefficient  structure  from  field  experiments,  even  in  the  simplest  of 
geological  setting  is,  therefore,  a  difficult  task. 

In  this  chapter  we  will  review  some  of  the  methods  proposed  in  the  literature 
for  obtaining  the  attenuation  profile  for  the  bottom  from  field  measurements.  We 
then  describe  a  field  experiment  that  has  been  proposed  by  Frisk  etal[6]  to  obtain 
the  plane-wave  reflection  coefficient  for  the  ocean  bottom  as  a  function  of  the  angle 
of  incidence.  The  sequence  of  processing  steps  that  can  be  adopted  to  obtain  the 
acoustical  properties  of  the  bottom  starting  with  the  plane-wave  reflection 
coefficient  information  is  then  indicated.  The  advantage  of  this  method  for  the 


Figure  3-1:  Configuration  of  experiment  by  Anderson  and  Blackman [7] 


o(f)  =  Attenuation  coefficient(dB/m)  of  sediment 

ow(f)  =  Attenuation  coefficient  of  water 

BL  as  Bottom  loss  of  bottom  reflected  signals  in  decibels 

D  =  Depth  of  water 

z  =  Thickness  of  sediment  wedge 

RLgb(f)  =  Peak  amplitude  of  sub-bottom  reflected  signal 

S(f)  =  Source  level  in  decibels 

SBL  =  Sub-bottom  loss 

TLg  —  Transmission  loss  of  signal  going  through  the  water /sediment 
interface 

TLw  as  Transmission  loss  of  signal  going  through  the  sediment  water 
interface 


determination  of  the  attenuation  profile  over  the  others  described  in  this  chapter 
will  also  be  presented.  We  conclude  the  chapter  by  a  statement  of  the  problem  that 
will  be  considered  in  the  remaining  part  of  this  thesis  and  the  assumptions  that  will 
be  made  in  solving  the  problem. 

3.1  Review  of  methods  in  the  literature 

Anderson  and  Blackman[7]  describe  a  method  for  determining  the  average 
attenuation  coefficient  from  insitu  measurements.  A  requirement  for  the 
experiment  is  that  the  sediment  layer  be  wedge-shaped  and  overly  a  highly 
reflective  bottom.  The  schematic  diagram  in  Figure  3-1  shows  such  a  configuration. 
In  such  an  environment,  the  field  experiment  is  performed  to  obtain  an  average 
value  of  attenuation  for  the  sediment.  Assuming  near  normal  incidence,  the 
pressure  pulse  due  to  the  source  interacts  with  the  water-sediment  interface  first. 
Part  of  the  energy  is  returned  to  the  water  column  while  the  remaining  is 
transmitted  into  the  sediment.  This  transmitted  energy  again  undergoes  reflection 
and  transmission  at  the  sediment  sub-bottom  interface.  Part  of  this  reflected 
energy  is  returned  to  the  water  column.  If  RLg  and  RLgb  are  the  received  pressure 
levels  of  the  sediment  and  sub-bottom  reflected  signals,  they  can  be  related  to  the 
source  level,  spreading  losses  and  other  losses  by  the  following  equations.  The 
difference  in  the  field  returned  from  the  water-sediment  interface  and  sediment- 
sub-bottom  interface  can  then  be  related  to  the  attenuation  in  the  sediment.  The 
symbols  used  in  the  following  equations  are  explained  in  the  legend  to  figure  3-1. 

Rlb  =  S-BL-20log2D-2awD  (3.1) 


RLgb  =  S-TLb-SBL-TLw-2az-20log(D+z)-2awD 


(3.2) 


Subtracting  ,  we  obtain 


RLb  -  RLsb  =  SBL  +  [TL^TLrBl^VHofiD\ 

+  2  az  +  20  log(D+z).  (3.3) 

Assuming  that  the  term  within  the  bracket  remains  constant  for  the  entire  wedge, 
we  obtain 

RLb—RLgb—20log2(D+z)—SBL  =  a(2  z)+c  (3.4) 

Ignoring  any  contribution  due  to  spreading  term,  i.e.  20  log2(D+z),  and  any  change 
in  the  sub-bottom  reflectivity,  this  equation  reduces  to 

RLb-RLeb  =  a(2  z)+c  (3.5) 

By  measuring  (RLb~RLgb)  at  different  points  in  range,  i.e.  at  different  depths 
in  the  bottom,  this  quantity  is  plotted  as  a  function  of  depth.  The  slope  of  the  line 
obtained  gives  the  attenuation  coefficient.  The  depth  at  each  observation  point  is 
determined  from  measuring  the  arrival  times  of  the  two  reflected  signals  and 
knowing  the  sound  speed  in  the  bottom. 

The  spectral  ratio  method  of  Jacobson,  Shor,  and  Dorman[8]  uses  seismic 
refraction  data  to  determine  attenuation  as  a  function  of  depth.  The  experimental 
configuration  is  shown  in  Figure  3-2.  Different  ray  paths  from  the  source  sample 
different  depths  in  the  ocean  bottom.  Assuming  horizontal  stratification,  the 
amplitude  spectrum  of  the  received  signal  waveform  can  be  written  as 


AJM*)  =  EQj(O,f)FJ(0J,Zm 


where  0  is  the  angle  of  incidence,/  is  the  frequency,  and  z  the  depth.A(0,/,z)  is  the 
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received  signal ,E(6,J)  is  the  source  signal,  f\0J,z)  the  earth  response  and  I(J)  the 
instrument  response.  The  suffix  j  refers  to  the  jth  observation.  The  earth  response 
function  can  be  written  as 

Fje,f,z)  =  TJkt,z)[:Jie,z)Gp,z)  (3.7) 

where  T  is  the  function  containing  transmission  and  reflection  coefficients, 
scattering  loss,  shear  conversion  etc.,  R  is  the  function  representing  spreading  loss 
and  G  is  the  attenuation  operator.  T  is  assumed  to  be  independent  of  frequency. 
The  frequency  dependence  of  G  is  given  by 

G(f,z)  =  exp(—k(z)JnPL(z)).  (3.8) 

where  PL(z )  is  the  path  length,  fk(z)  is  the  attenuation  coefficient.  If  E  and  /  are 
known,  the  earth  response  F  can  be  obtained.  To  get  a  good  estimate  of  the  source 
function,  the  direct  water  borne  energy  is  measured.  This  is  given  by 


=  E0m J0,zwm 


where  /?  •  is  the  spreading  loss  function  ralated  to  the  waterborne  path.  The 
spectral  ratio  defined  as  the  logarithm  of  the  ratio  of  the  refracted  signal  amplitude 
to  the  direct  water  signal  is  then  obtained. 


A(0,/,*)  Tp,t)R]fitzY2[r%») 

Taking  the  logarithm  of  both  sides,  the  spectral  ratio  is 


(3.10) 


SR(e,f,z)  =  lnTj(9,z)+lnRj(0,z)— 

InR^+lnG^z) 


(3.11) 


Assuming  that  the  attenuation  cofficent  is  linearly  related  to  frequency 


lnG(f,z)  =  -k(z)fPL{z) 


(3.12) 


Then, 


SRjDJ.z)  =  SRfi ,(«,*)-£  k%JPLtj 


(3.13) 


where, 


SRfi  =  lnTj,0,z)+lnRj,0,z)-lnRwj(0,z) 


(3.14) 


By  plotting  SR(9,f,z )  with  respect  to  /  we  can  obtain  the  intercept  SRQ.  Then, 


SRr-SR 


■f—*-  -Ev% 


(3.15) 


where  SRQ  is  the  intercept,  PL-  is  the  path  length  in  the  ith  layer  for  the  ;th 
observation.  A  linear  inverse  method  is  now  used  to  obtain  Having  obtained  k ., 
the  attenuation  in  the  ith  layer  equals  fk-.  Wrolstad[9]  also  uses  a  similar  method 
for  determining  the  attenuation  profile.  The  spectral  ratio  method  has  been  applied 
by  Hauge[10]  to  vertical  seismic  profiles  to  obtain  attenuation  information. 

Stoll  and  Houtz(ll]  utilize  the  fact  that  in  certain  areas  of  the  ocean  where 
there  are  deep  sediment  layers,  the  velocity  profile  is  smooth  and  linear.  With 
linear  profiles  the  refracting  paths  in  the  sediment  layers  will  be  arcs  of  a  circle. 
The  experiment  is  similar  to  a  refraction  experiment  and  the  configuration  is  shown 
in  Figure  3-3.  The  amplitude  of  the  signal  received  at  the  time  of  the  first  arrival  or 


at  a  multiple  of  this  time  is  given  by 


1  tan(0)  . 

\=\w  — (— (3.i6| 
l+n  S(x,&) 

where, 

dx 

S(x,0)  =  x—  (3.17) 

dO 

where  n=0  for  the  first  the  first  arrival  and  n=l  ,2,..  for  the  first  multiple  and  so 
on.  The  term  in  parenthesis  accounts  for  geometrical  spreading, R(6)  is  the 
reflection  coefficient  and  .A(0)  is  the  amplitude  of  the  source.  9  is  the  angle  of 
inclination  of  the  angle  at  the  source  and  (l-fn)x  is  the  range  between  source  and 
receiver.  The  exponent  -  d>(6)  is  the  total  loss  of  amplitude  occuring  due  to  intrinsic 
attenuation,  scattering  and  other  causes  and  may  be  expressed  in  terms  of  the 
attenuation  coefficient  by  the  expression  given  below. 

4>n  =  (n+l)J^  a(J)As  (3.18) 

S 

To  evaluate  0,  we  form  the  ratio  AJA^, 


/?  =  —  = 

Therefore, 


(3.19) 


2Ra 

■fn  =  — ln( — r)  (3.20| 

0  (i-«V 

and  knowing  $Q,  a(f)  is  determined. 

Mitchell  and  Focke[12]  use  the  setup  shown  in  Figure  3-4  for  their  field 


experiment.  One  of  the  ray  paths  from  the  source  to  the  receiver  is  shown  in  the 
Figure  3-4.  There  may  be  many  more  such  paths.  The  major  part  of  the  loss  of 
signal  energy  occurs  along  portions  of  the  ray  that  lie  in  the  sediment.  This  loss, 
when  normalised  to  a  single  bottom  encounter  and  expressed  as  a  function  of 
grazing  angle  is  the  bottom  reflection  loss.  Obtaining  this  is  the  first  step  in  this 
method.  Mitchell  and  Focke  then  use  the  fact  that  for  rays  penetrating  the  bottom 
at  high  angles  of  incidence,  the  bottom  loss  associated  with  them  is  largely  due  to 
attenuation  in  the  layers.  Under  such  assumptions,  the  bottom  loss  can  be 
expressed  as  a  function  of  the  attenuation  coefficient  and  the  path  length  in  the 
bottom.  When  the  attenuation  is  a  function  of  depth,  then  this  relationship  is 
given  by  the  equation  below. 


=  rhm 


Here  we  have  assumed  that  a(f,z)=k(z)f  and  the  integral  is  over  the  path  length. 
Using  this  and  considering  a  number  of  rays  all  at  large  angles  of  incidence  a  matrix 
of  the  form  given  below  can  be  constructed. 


(3.22) 


Here  L is  the  path  length  in  the  jth  layer  of  a  ray  at  an  angle  of  incidence  of 
•.  To  compute  the  path  length  it  is  assumed  that  the  sound  speed  in  the  bottom  is 
known.  This  matrix,  which  is  lower  triangular,  can  be  easily  inverted  to  give  k .  in 
different  layers  and  as  before  the  attenuation  coefficient  in  each  layer  can  thus  be 
obtained. 


Spofford[13]  describes  a  method  of  determining  the  sound  speed  gradient  and 
an  average  value  of  the  attenuation  coefficient  for  the  sediment  from  a  knowledge 
of  the  bottom  loss  expressed  as  a  function  of  the  angle  of  incidence.  Such  data 
exhibit  an  abrupt  increase  in  loss  at  a  grazing  angle  corresponding  to  the 
development  of  a  minimum  range  caustic  in  the  bottom  refracted  path.  If  this 
grazing  angle  is  obtained  from  the  data,  then  using  a  relationship  between  this 
angle  and  the  sound  speed  gradient  for  the  bottom,  the  sound  speed  gradient  is 
determined.  The  knowledge  of  the  sound  speed  gradient  gives  the  path  lengths  of 
the  refracted  rays.  The  attenuation  coefficient  value  is  then  obtained  from  the 
difference  in  levels  between  the  bottom  reflected  and  refracted  rays. 

The  methods  described  above  make  some  or  all  of  the  following  assumptions. 

1.  The  reflection  and  transmission  coefficients  are  independent  of 
frequency. 

2.  The  attenuation  coefficient  is  linearly  related  to  frequency  and 
therefore  a(f,z)  can  be  written  as  a(f,z)—k(z)f. 

3.  The  compressional  wave  speed  for  the  medium  is  known. 

4.  The  compressional  waves  are  dispersionless. 

5.  Where  other  loss  mechanisms  like  scatter, shear  conversion  are 
involved,  the  attenuation  coefficient  obtained  will  be  the  ’effective 
attenuation’  and  not  the  'intrinsic  attenuation’. 

The  field  measurement  of  attenuation  coefficent  show  considerable  scatter  as  seen 
in  Figure  2-11.  On  the  basis  of  the  data  presented  in  this  figure, it  is  clear  that 
extrapolation  of  the  high  frequency  results  into  the  low  frequency  region  using 
linear  law  will  be  in  error.  Further,  the  data  in  the  low  frequency  range  are  spread 
over  a  wide  range  of  values  of  attenuation  and  there  is  no  consistent  pattern  to 
justify  a  particular  frquency  dependence. 


(4.8) 


dz 


Where  k2  is  the  wave  number  in  region  HI  (i.e.)  fc2=aj/C2. 

For  an  incident  wave  of  unit  amplitude,  the  field  in  region  I  is  the  sum  of  the 
incident  and  reflected  waves. 

P0=etkzOz+R(kx)e~ikzOz  (4.9) 

In  the  above  equation  R(kx)  is  the  plane  wave  reflection  coefficient  of  the  entire 
bottom  structure  and  kzQ=  (fr02  —  fcx2)^2  is  the  vertical  wave  number  in  the  water 
column.  The  subbottom  is  an  infinite  half  space  and  only  outward  propagating 
wave  exists  in  this  region.  The  field  in  the  subbottom  is,  therefore,  given  by  the 
following  expression. 

P2=T\kx)eikz2z  (4.10) 

where  kz2  is  the  vertical  wave  number  in  the  subbottom  and  T \kx)  is  the 
transmission  coefficient.  Note  that  kz2  is  complex  to  allow  for  attenuation  in  the 
bottom.  The  pressure  field  in  the  sediment  layer  is  obtained  by  solving  equation 
(4.7).  The  boundary  conditions  that  the  solution  has  to  satisfy  are  the  continuity  of 
pressure  and  normal  particle  velocity  at  the  two  interfaces  i.e  at  2=0  and  z=h 
representing  the  water-sediment  and  sediment-subbottom  interfaces.  The 
application  of  these  conditions  at  the  sediment-water  interface  leads  to  the 
following  equations. 


fl(0)=l+/? 


(4.11) 


(4.2) 


dz 2 


M*)( 


1  dP[z) 

—)' - 

p(z)  dz 


(k-(z)-k*)P[z)=0. 


Defining  a  new  variable  v(z)  as  v(z)=p  '^(aj/^a)  and  substituting  in  the 
equation  we  obtain  the  following  relation. 


d2v(z) 

dz*~ 


±Q(z)v(z)=0 


where, 


p'l'\z)ft(z\ 

2  yi$ 


When  the  density  in  the  medium  is  a  constant  the  wave  equation  takes  the 
given  below. 


d2P[z)  n 

— — +(t2(a)-*I2)fla)=0 
dz 

The  wave  equation  in  the  three  regions  will  then  be  as  given  below. 
Water  column  (Region  1): 


d-P, 


2+(*o2-*IVo= 0 


dz* 

where  kQ  is  the  wave  number  in  region  I  (i.e.)  k0=u/CQ. 
Sediment  layers  (Region  II): 


d2 

—o+Q(z)v(z)=0 
dz “ 

Sub-bottom  (Region  III): 


*•  k’w  , 


wave 


(4.3) 


(4.4) 

form 


(4.5) 


(4.6) 


(4.7) 
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Figure  4-1:  Ocean  bottom  model 


Chapter  4 

The  Forward  Problem 


The  forward  problem  deals  with  the  determination  of  the  response  of  a 
system  to  a  known  input  given  system  parameters.  In  the  context  of  our  problem,  it 
deals  with  the  determination  of  the  acoustic  field  in  the  sediment  layers  for  a  plane 
wave  of  unit  amplitude  incident  on  the  top  of  the  sediment.  The  acoustic 
characteristics  of  the  sediment  layer  and  the  sub-bottom  are  assumed  known.  We 
need  to  solve  the  forward  problem  for  two  reasons.  First,  the  plane  wave  reflection 
coefficient  and  the  pressure  field  in  the  sediment  layers  for  the  assumed  guess 
model  are  needed  to  solve  the  inverse  problem  as  we  will  see  in  the  next  chapter. 
Second,  we  use  synthetic  data  to  test  the  inversion  scheme  and  for  this  purpose  the 
plane  wave  reflection  coefficient  is  required  to  be  computed.  Two  methods  for 
solving  the  forward  problem  are  presented  in  this  chapter;  a  uniform  asymptotic 
method  and  a  method  based  on  a  propagator  matrix  formulation. 

We  use  the  ocean  bottom  model  developed  in  Chapter  2  and  represented  in 
Figure  4-1.  Consider  a  plane  wave  of  unit  amplitude  incident  at  the  water  sediment 
interface  inclined  at  an  angle  6  to  the  vertical.  The  acoustic  field  in  the  ocean  and 
the  bottom  obey  the  acoustic  wave  equation.  If  <f>  represents  the  pressure  field;  the 
wave  equation  for  <p  is 


d2<f>  d2<j>  1  d  <j>  n 

— +— +/>W(— i'— +fc2<2>*=°- 
ox*  o  z*  p{z)  a  z 


(4.1) 


Writing  Mx,z)—P[z)eikzx  the  wave  equation  becomes 
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circular  symmetry  exists  m  the  point  source  field. 

3.  The  sediments  can  be  modelled  as  a  fluid. 

To  simplify  the  analysis,  we  further  assume  that  the  acoustic  parameters  are 
known  in  the  region  z  <  0  and  z  >  h  and  that  the  parameters  need  be  determined 
only  in  the  region  0  <  z  <  h. 

/e  start  by  assuming  that  there  are  no  errors  in  the  sound  speed  and  density 
and  on! y  the  attenuation  profile  needs  to  be  determined.  Then,  treating  errors  in 
the  sound  speed  and  density  as  unknowns,  the  problem  is  extended  to  the 
determination  of  these  as  well. 
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Figure  3-7:  Exact  and  guess  model  for  the  ocean  bottom 


1.  Since  a  C.W.  source  is  used  the  effect  of  dispersion  is  not  relevant.  The 
sound  speed  obtained  in  step  2  of  the  scheme  is  the  value  at  the 
frequency  of  the  source. 

2.  No  a  priori  assumption  is  made  on  the  frequency  dependence  of 
attenuation.  The  scheme  provides  the  attenuation  as  a  function  of 
depth  at  the  operating  frequency 

3.  By  doing  simultaneous  inversion  for  attenuation  coefficient  and 
corrections  to  sound  speed  and  density  structure  ,  the  effects  of  errors 
in  sound  speed  and  density  on  the  attenuation  profile  can  be  overcome. 


3.3  Statement  of  the  problem  for  the  determination  of  attenuation 
coefficient 

For  the  purpose  of  this  dissertation  we  will  assume  that  by  some  suitable 
experiment  and  subsequent  processing  the  plane  wave  reflection  coefficient  as  a 
function  of  angle  of  incidence  has  been  obtained.  Further,  we  assume  that  this 
information  has  been  used  in  an  inverse  scheme  to  obtain  a  close  estimate  of  both 
compressional  wave  speed  and  the  density  profile  for  the  bottom. 

Consider  the  two  models  of  the  ocean  bottom  in  Figure  3-7.  The  model  on  the 
top  is  the  true  model  of  the  ocean  bottom  with  density,  sound  speed  and 
attenuation  profiles  as  indicated.  The  other  model  is  our  present  knowledge  of  the 
parameters  describing  the  ocean  bottom(i.e)  a  close  estimate  of  sound  speed  and 
density  stucture  as  obtained  in  the  earlier  stage  of  processing.  The  problem  is  to 
develop  a  method  to  determine  the  attenuation  profile  and  corrections  to  the  sound 
speed  structure  and  density  structure  for  the  sediment  layer.  In  solving  the  problem 
the  following  assumptions  will  be  made. 

1.  The  ocean  bottom  is  horizontally  stratified. 


the  pressure  field  by  the  foilwing  equation. 

—R{kr){exp{ikzQ{zr+zg))}  =  f  PreJ(r)J0(krr)rdr  (3.24) 

kz0  Jo 

Therefore,  by  performing  the  Hankel  transform  of  the  reflected  pressure  field 
we  can  obtain  the  plane  wave  reflection  coefficients  as  a  function  of  the  horizontal 
wave  number  or  equivalently  as  a  function  of  the  angle  of  incidence.  Suitable 
methods  have  been  developed  for  performing  the  Hankel  transform 
operation[14-18].  Schoeberg(19]  uses  a  different  approach  to  obtain  the  plane-wave 
reflection  coefficient  from  a  measurement  of  the  pressure  field. 

3.2.2  Determination  of  sound  speed  and  density  profiIes(step  2) 

In  this  step  the  sound  speed  and  density  profiles  are  obtained  using  the  plane- 
wave  reflection  coefficients  as  input  data.  Direct  inverse  methods [20-22]  can  be  used 
for  this  purpose. 

3.2.3  Determination  of  attenuation  coefllcient(step  3) 

Using  the  information  on  sound  speed  and  density  values  obtained  in  the  last 
step  and  the  information  on  the  plane  wave  reflection  coefficient,  the  attenuation 
coefficient  profile  is  obtained.  In  this  dissertation  a  perturbation  approach  is  used 
to  accomplish  this. 

This  method  of  determining  the  attenuation  coefficient  profile  from  plane- 
wave  reflection  coefficient  obtained  with  a  CW  source  experiment  has  the  following 
advantages  over  the  other  methods  where  a  broadband  source  is  used. 
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pressure  field  by  studying  the  peak  to  null  difference  m  the  measured  interference 
pattern. 


3.2  Determination  of  the  attenuation  profile  from  the  plane-wave 
reflection  coefficient 


By  suitably  proceessing  the  measured  pressure  field,  the  plane-wave  reflection 
coefficient  for  the  ocean  bottom  can  be  obtained.  The  plane-wave  reflection 
obtained  is  a  function  of  angle  of  incidence  and  it  contains  all  the  information 
about  the  acoustic  parameters  of  the  bottom.  The  determination  of  the  acoustical 
properties  of  the  bottom  from  the  plane-wave  reflection  coefficients  is  a  more  exact 
method  and  we  will  follow  this  path. 

A  possible  sequence  of  steps  to  obtain  the  complete  set  of  acoustic  parameters 
starting  with  the  measured  field  is  given  in  Figure  3-6.  We  will  briefly  indicate  the 
approaches  that  can  be  adopted  for  solving  each  step. 


3.2.1  Determination  of  plane-wave  reflection  coefficient(step  1) 


The  total  pressure  field  is  the  sum  of  the  direct  field  and  the  reflected  field. 
After  subtracting  out  the  direct  field,  the  reflected  field  is  given  by[14]; 


PreJ(r)  =  f  —R( kr)  JQ{  kr)  { exp{  ikj  ^^+  zg))}  kd  kf 
Jo  k.n 


(3.23) 


where  R(kf)  is  the  plane  wave  reflection  coefficient  for  the  bottom,  kf is  the 
horizontal  wave  number,  r  is  the  range,  kgQ is  the  vertical  wave  number  in  the  wate  r 
column  and  zr  and  zg  are  receiver  and  source  heights  from  the  bottom.  JQ  is  the 
zeroeth  order  Bessel  fuction.  Alternatively  we  can  relate  the  reflection  coefficient  to 
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Figure  3-5:  Configuration  of  experiment  by  Frisk  etal[6] 


We  noted  in  Chapter  2  that  both  the  theoretical  models  predict  substantial 
dispersion  for  sediments  with  high  attenuation  coefficent.  In  view  of  this  the 
assumption  that  the  medium  is  dispersionless  is  questionable.  This  is  an  important 
consideration  because  all  methods  described  above  assume  that  the  compressional 
wave  velocity  profile  for  the  medium  is  known  and  this  information  is  used  in  the 
determination  of  the  attenuation  coefficient.  Therefore  any  errors  in  the 
compressional  wave  speed  arising  out  of  this  assumption  will  manifest  as  errors  in 
attenuation  coefficient. 

These  considrations  indicate  the  desirability  of  an  experiment  which  uses  a 
monochromatic  source  for  the  determination  of  the  acoustic  parameters  of  marine 
sediments.  Experimental  configuration  of  one  such  field  measurement  carried  out 
by  Frisk  etal  is  shown  in  Figure  3-5.  Two  receivers  are  moored  on  the  ocean  bottom 
at  heights  1.17  m  and  54.55  m  from  the  bottom.  A  CW  source  is  suspended  from 
the  ship  at  the  end  of  a  long  cable.  The  ship  is  then  allowed  to  drift  and  as  the  ship 
drifts  away  slowly  the  two  receivers  record  the  pressure  field.  By  suitable  signal 
processing  methods  the  pressure  field  as  a  function  of  range  is  obtained. 

Since  the  recorded  pressure  field  at  the  receiver  depends  on  the  acoustic 
parameters  of  the  bottom,  these  parameters  can  be  obtained  by  forward  modelling. 
In  this  we  assume  an  initial  guess  model  for  the  acoustic  parameters  and  solve  the 
direct  problem  to  obtain  the  pressure  field  at  the  receiver  as  a  function  of  range. 
This  is  compared  with  experimentally  determined  value,  and  the  acoustic 
parameters  are  changed  until  good  fit  between  the  experimentally  determined  value 
and  the  theoretically  determined  value  is  obtained.  Frisk[6]  has  done  such  analysis 
with  considerable  success  and  obtained  geo-acoustic  models,  including  an  average 
attenuation  coefficient,  for  different  regions  of  the  ocean.  Frisk[6)  also  describes 
how  an  estimate  of  the  attenuation  coefficient  is  obtained  from  the  recorded 


- P'(0)= — —{1— i?) 

P(0)  P0 


(4.12) 


In  the  above  equations  p(0)  is  the  density  in  the  sediment  at  the  interface  while  p0 
is  the  density  in  the  water  column.  P(0)  is  the  pressure  field  in  the  sediment  at 
2=0.  The  primes  denote  derivatives  with  respect  to  z.  Similarly  applying  the 
boundary  conditions  at  the  other  interface  yields  a  second  set  of  equations. 


P[h)=Teikz2h 


-P  '(h)=—^—-€*kz2h 


(4.13) 


(4-14) 


4.1  Uniform  asymptotic  solution 

The  uniform  asymptotic  solution[l]  to  the  wave  equation  in  region  II  has  been 
discussed  in  detail  by  Kawahara[2].  The  solution  that  is  valid  when  the  function 


Q(z)  varies  smoothly  and  there  is  only  one  turning  point,  is 


u(z)=c1Q-I/4(zh,/4Aj(T(Z)H-c2Q-|/4(zh'/<Bj(T(2)) 
where  A j  and  Bj  are  the  Airy  functions.  t(z)  is  given  by  the  expression 


(4.15) 


iw=3/2(  f V/2mi2/3 


(4.16) 


and  zt  corresponds  to  the  turning  point  i.e  Q(zt)=0.  The  solution  for 
P[z )  =  pl!2{z)v(z)  is  then; 


Pz)=cigi{z)+c2g2(z) 


(4.17) 


j2(z)=C2//2(z)e-,/4(z)'r1/4B1«z)) 

and  Cj  and  c2  are  arbitrary  constants. 

Applying  the  boundary  conditions  at  the  two  interfaces  we  now  obtain  a  set 
of  four  equations  with  four  unknowns  to  be  determined  namely  the  cofficients  Cj 
and  c2  and  the  reflection  and  transmission  coefficients  R  and  T.  Solving  for  these 
unknowns,  we  obtain  R^T^^andc^-  The  coefficients  are  now  substituted  in  equation 
(4.17)  and  the  field  obtained.  Kawahara[2]  has  shown  that  for  certain  canonical 
profiles  like  c  linear,  c2  linear  and  c3  linear  the  integral  in  (4.16)  can  be  performed 
analytically. 

4.2  Propagator  matrix  method 

The  second  method  for  solving  the  forward  problem  models  the  sediment  as  a 
stack  of  homogeneous  layers.  By  making  the  thickness  of  the  layers  small  any 
arbitrary  variation  of  density,  sound  speed  and  attenuation  can  be  accommodated. 
Mook[3)  has  proposed  a  method  for  determining  the  reflection  coefficient  using  the 
propagator  martix  approach.  We  will  essentially  follow  his  approach  and  indicate 
how  a  numerically  stable  algorithm  is  obtained  for  computing  the  plane  wave 
reflection  coefficients  and  the  field  in  the  sediment  layer. 

The  propagator  matrix  method  is  dealt  with  extensively  in  the  literature[4-6|. 
We  will  explain  in  brief  the  essential  features  of  the  propagator  matrix  method 
before  taking  up  the  problem  indicated  above. 


Any  second  order  differential  equation  can  be  written  as  a  set  of  first  order 
equations.  The  general  representation  is  given  below. 


df 

- A(z)f(z)=0 

dz 


(4.18) 


For  the  equation  (4.2),  f  is  an  2x1  vector  and  A  is  a  2x2  matrix.  The  propagator 
matrix  (also  called  matrizant)  is  defined  as 


p(z,z0)  =  i+  f A(e1)de,+ 

J*0 

•Mn)  JZn 


(4.19) 


where  I  is  the  identity  matrix  of  the  order  of  A.  By  substitution  it  can  be  shown 
that  P(z,zQ)  satisfies  the  following  differential  equation 


d2P{z,zQ) 


-  A(z)P(z,z0)=0 


(4.20) 


Also  P(z0,Zq)=I.  The  most  important  property  of  the  propagator  matrix  is  that  it 
enables  us  to  obtain  f(z)  from  its  value  at  some  other  point  Zq. 


t[z)  =  P(z,z0)f(z0) 


To  see  that  this  is  true, we  substitute  f(z)  in  the  differential  equation.  Then, 


(4.21) 


d2i[z) 


—  A(z)f(z)=0 


(4.22) 


Therefore,  f(z)  satisfies  the  differential  equation.  Also  it  gives  the  value  of  f(z) 


at  zQ  correctly. 
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f(20)=P(^0>2o)f(2o)=f(Zo) 


(4.23) 


Another  interesting  property  of  the  propagator  matrix  is  shown  by  the  following 
equations. 


f(’2)  =  p(z2>2i)f(2i) 

=P(z2,zl)P(zvzQ)f{z0) 


(4.24) 


If  we  now  choose  z„  =  zQ  then, 


f(-’0)  =  P(Z(].2,)P(-1i?0lf<J0l 


(4.25) 


I  =  P(20'zl*P(2l’zol 


(1.26) 


When  the  matrix  A  is  not  dependent  on  z  as  in  the  case  of  homogeneous 


layers, 


P(^,2o)=i+(2-^o)a+1/2(^-20)2a2+ 


(4.27) 


P(2,0o)=e("_2o)A 


(4.28) 


For  a  square  matrix  A  with  distinct  eigenvalues,  \k(k=l,  •  •  •  ,n)  by  the  Sylvester 
theorem, 


n a)  =  E  II 


,  ^k~^r 


(4.20) 


jfc=l  rjLk  *  T 

Reverting  now  to  the  differential  equation  (4.2)  it  is  written  as  a  set  of  first 


order  equations. 


When  applied  to  a  stack  of  homogeneous  layers  p( z)  and  kz  represent  the  density 
and  the  vertical  wave  number  in  the  layer  under  consideration.  As  before  u 
represents  the  frequency  and  the  vector  f  contains  the  elements  P  and  w,the 
pressure  and  normal  particle  velocity.  The  matrix  A  for  the  above  equation  is, 
therefore 


=  **  \k 


up(z) 


up(z) 


The  propagator  matrix  is  now  obtained  using  equation  (4.29). 


(4.31) 


P(»*0> 


cos  (kJz-Zg)) 
i 

— k/mfz-z^ 
up 


\ up 

—sin  (z-Zq) 

kz 

cos(z-Zq) 


(4.32) 


In  the  above,  Y==—  is  the  admittance. 
up 

Consider  a  layer  of  thickness  h  as  shown  in  figure  4-2.  The  pressure  and 
particle  velocity  at  depth  z= 0  is  related  to  the  pressure  and  particle  velocity  at 
depth  z—h  by  the  following  equation. 


c°s  (kzJh) 


Y(l) 

iY[l)sin(kzth)  cos  (ktJh) 


\n(k  h) 


(4.33) 


Since  each  layer  is  homogeneous,  the  pressure  and  normal  particle  velocity  in 
the  layers  can  be  written  as  a  sum  of  up  and  down  going  waves. 


P[z)  =  Deikzl\2~zrejl  +  Ue~ikz\\z-zrefl 


(4.34) 


w{z)  =  —(De'kzl\z~zre/  -  Ue~ikzl\z-zre/) 
up 


(4.35) 


P(z)  1  1 

—  fr  mjr 

w(2)\  [—  ~ 

up  up 


elkzllzzrejl  0 
0  eikzlfzzre j 


Setting  zrej=  0,  we  obtain 


(4.36) 


eikzlz 


^l)e'V 


eikzlz 


,-i  k_,z 


Vllje-1^/ 


(4.37) 


where, 


kz\ 

Hi)=— 


Now  consider  a  stack  of  homogeneous  layers  as  shown  in  Figure  4-3.  Then  the 
pressure  and  particle  velocity  at  z=0  is  related  to  the  up  and  down  going  waves  in 
the  watercolumn  by  the  following  equation. 


Y(0)  ->10)  111/(0) 


(4.38) 


or  equivalently, 


Following  the  method  of  Mook[3]  the  matrix  in  equation  (4.33)  is  modified 
and  the  elements  expressed  in  terms  of  ratio  of  impedances.  This  is  done  to  avoid 
difficulties  that  may  arise  due  to  P  and  w  being  of  different  scales  of  magnitude. 


rfoMO) 


cos  (kjl)h)  -bin  (k2(l)h)  P(l) 

-i£(l)sin(kjl)h)  £{l)cos  (kjl)h) 


(4.40) 


where  £(l)=y(l)/Y(0),  f(0)=l/Y{0),  and  ff(z)(l)=l/y(l).  Using  the  propagator 
matrix  we  can  relate  the  P  and  w  at  the  nth  interface  to  the  values  of  P  and  w  at 
the  first  interface: 


rfOMO) 


A  n-D 

c(n-lMn-l) 


(4.41) 


where  <f> s  are  the  propagator  matrices  for  the  layers  from  1  to  n-1.  The  pressure 
and  particle  velocity  at  the  bottommost  interface  is  now  related  to  the  up  and 
down  going  waves  in  the  terminating  half  space. 


P(n-1) 

c(n-l)u<n-l) 


1  1 


fln)  -f(n)J[U(n) 


(4.42) 


The  relation  between  D{0),U(0)  and  D(n),U(n)  can  now  be  obtained. 


AO)  i  1  1  11  An) 

=  -  (*, . «„.!>  H«) 

MO)  2  1  -1  «■>)-«(»)  l\  n) 


=  <p 


D(n)+U[n) 

HnMnyU{n)} 


(4.44) 


where 


1  1  1 

<P  =  - 

2  1  -1 


(^1 . 


(4.45) 


and  <pis  a  2x2  matrix  with  elements  0n'^i2’^2lan<^22‘  In  the  infinite  half  space 
representing  the  subbottom  U(n)=0.  Therefore, 


(4.46) 


m 

=  <t> 

D{  n) 

m 

The  reflection  coefficient  is  then  given  by 

M»)  ^2i+fln^22 

R  =  -  =  - 


(4.47) 


D(0)  ^n+£(«)^12 

Having  obtained  the  reflection  coefficient  the  pressure  and  particle  velocity  at 
the  top  interface  are  easily  found  to  be  fl(0)=l+i?  and  u^0)==  V(0)(1— /?) 
respectively.  We  can  now  propagate  the  pressure  and  particle  velocity  downward 
and  determine  it  at  any  layer.  This  is  similar  to  solving  the  differential  equation 
given  its  initial  conditions.  If  this  procedure  is  adopted,  the  scheme  becomes 
numerically  unstable  ami  the  solution  rapidly  diverges  from  the  correct  solution  in 
regions  where  the  waves  become  inhomogeneous.  Hawker  and  Foreman[7]  give 
examples  to  show  that  the  direction  in  which  the  integration  is  carried  out  is 
important  to  obtain  a  stable  solution  and  recommend  that  the  integration  be  done 


upwards  starting  with  the  conditions  at  the  bottommost  interface  as  the  initial 
conditions. 

To  do  this,  we  assume  that  Z)(n)=1.0  and  U(n)=0.  The  pressure  and  particle 
velocity  at  the  nth  interface  is  then  related  to  D(n)  and  U(n). 


n  n-l) 

{(n-l)t^n) 


1  1  If1 

^(n)  -£(n)  0 


(4.48) 


We  now  relate  /’(n-l)  and  u(n-l)  to  P(n-2)  and  u(n-2). 

P[n-2)  1  \f\n-l) 

=  <f>  (4.49) 

f(n*2)u(n-2)  f(n-l)u>(n-l) 


As  before  relating  the  up  and  down  going  waves  at  the  top  to  the  up  and  down 
going  waves  in  the  infinite  halfspace  in  the  bottom,  we  obtain  the  following 
expression. 

m  [i 

=  $  (4.50) 

.MO)  U(n)  . 


This  gives  the  reflection  coefficient  as 


_  MW  _  ^21+£W22 

MO)  <t>n+t{n)<t>l2 

The  pressure  computed  will  correspond  to  an  incident  amplitude  of  MO)  The 
pressure  field  at  any  point  in  the  layers  is  then  obtained  by  dividing  thro  gh  by 
MO).  In  this  thesis  we  use  the  propagator  matrix  method  to  compute  the  reflection 
coefficient  and  the  pressure  field  in  the  sediment  layers. 
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Chapter  5 

Formulation  of  the  Inverse  Problem 


In  Chapter  3  we  proposed  a  sequence  of  steps  for  the  determination  of  the 
acoustic  parameters  from  the  measured  pressure  field.  We  also  gave  a  statement  of 
the  problem  that  is  addressed  in  this  thesis.  In  Chapter  4  we  dealt  with  the  solution 
of  the  forward  problem.  In  this  chapter  we  formulate  the  inverse  problem,  namely, 
the  determination  of  the  attenuation  coefficient  profile  in  the  sediment  layer  using 
the  plane  wave  reflection  coefficient  at  various  angles  of  incidence  as  input 
information.  The  real  and  imginary  part  of  the  wavenumber  at  220  Hz  are  given  in 
Table  5-1  for  different  types  of  marine  sediments.  These  have  been  obtained  using 
data  presented  by  Hamiltonjl]  for  the  compressional  wave  speed  and  attenuation 
coefficient.  The  attenuation  coefficient  at  220  Hz  is  computed  by  extrapolating  into 
the  low  frequency  region  using  a  linear  law.  The  attenuation  coefficient  is  order  of 
magnitude  smaller  than  the  real  part  of  the  wavenumber  and  can  therefore,  be 
treated  as  a  perturbation  of  the  real  part.  Hence,  we  formulate  the  inverse  problem 
largely  following  the  perturbation  method  introduced  by  Cohen  and  Bleistein[2|. 
They  dealt  with  the  problem  of  determining  perturbations  in  sound  speed  structure 
using  the  plane  wave  reflection  coefficient  at  all  frequencies  and  at  a  fixed  angle  of 
incidence  as  input  information.  In  Chapter  3  we  indicated  the  advantage  of  using  a 
constant  frequency  source  for  determining  the  acoustic  properties  of  marine 
sediments.  Therefore,  an  essential  difference  between  our  formulation  and  theirs  is 
that  we  use  plane  wave  reflection  coefficient  at  different  angles  of  incidence  and  at 
a  fixed  frequency  as  input  information. 
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In  this  chapter,  we  first  arrive  at  a  non-linear  integral  equation  relating  the 
plane  wave  reflection  coefficient  to  the  unknown  attenuation  coefficient  profile.  The 
integral  equation  is  then  linearised  using  the  Born  approximation.  We  also  obtain 
the  region  of  validity  of  the  Born  approximation. 

Rytov  approximation  can  also  be  used  to  obtain  a  similar  linear  integral 
equation.  Keller[3]  has  shown  that  the  Rytov  approximation  is  valid  for  longer 
ranges  than  the  Born  approximation  in  regions  where  there  is  only  one  field,  for 
example  in  regions  where  only  transmitted  field  exists.  Keller  observes  that  where 
there  are  more  than  one  field,  the  Rytov  method  has  to  be  applied  to  each  wave 
separately  for  it  to  retain  its  advantage  over  the  the  Born  approximation.  On  the 
other  hand,  the  Born  approximation  can  be  applied  to  the  total  field.  In  a  recent 
work  Oristaglio[4j  studied  the  accuracy  of  the  Born  and  Rytov  approximation  to 
the  laws  of  reflection  and  refraction  and  noted  that  the  Rytov  method  works 
reasonably  well  for  in  the  backscattering  regions  as  well  even  when  the  method  is 
applied  to  the  total  field;  a  surprising  result  in  the  context  of  Kellers  observations. 
We  formulate  the  inverse  problem  using  the  Rytov  method  so  as  to  make 
comparisons  between  the  two  methods. 

Consider  the  two  models  in  Figure  3-7.  One  of  these  is  the  true  model  and  the 
other  the  guess  model  which  we  call  the  background  model.  The  background  model 
represents  our  present  knowledge  of  the  model.  In  this  chapter  we  assume  that  the 
only  unknown  is  the  attenuation  coefficient  profile  in  the  sediment  layer.  The  other 
acoustic  parameters  of  the  two  models  are  the  same.  We  first  determine  the 
response  at  the  water  sediment  interface  of  the  two  models  to  a  plane  wave  of  unit 
amplitude.  For  the  true  model  this  is  obtained  from  a  field  experiment.  For  the 
background  model  we  obtain  this  by  solving  the  forward  problem  using  the  method 
described  in  Chapter  4.  We,  then,  relate  the  difference  in  response  to  the  unknown 


attenuation  coefficient  profile  and  obtain  an  integral  equation  representation  of  the 
problem. 


5.1  Inverse  problem  formulation  using  Born  approximation 

5.1.1  Derivation  of  integral  equation 


Consider  a  plane  wave  of  unit  amplitude  incident  at  an  angle  0  as  shown  in 
figure  5-1.  The  pressure  field  in  the  sediment  layer  is  given  by  the  equation; 


d2<f>  d2<f> 


d<f> 


a  2  ■  2+P( 2)(“)'t~+  k2(z)4>=0 

dx2  dz2  p(z)  d  z 


(5.1) 


Let  <(>(x,z)  —  P[z)exkxx.  Substituting  in  the  equation  (5.1)  we  obtain  the 
following  one  dimensional  wave  equation  for  the  field  in  the  sediment. 


d2P[z )  1  dP[z)  „ 

-5— +o(^(— c- — (5.2) 
dz£  p(z)  dz 

Let  v(z)  =  p~xl2(z)P[z).  Substituting  this  in  equation  (5.2),  we  obtain  the 
Schrodinger  type  equation  given  below. 


where, 


(5.3) 


p't-{z)  f(z)  , 
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We  derive  a  similar  equation  for  the  backgroud  model  as  well.  The  subscript 


’b’  in  the  following  equations  refer  to  the  backgroud  model. 
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Figure  5-1;  Ocean  bottom  model 
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The  wavenumber  in  the  sediment  layer  is  a  complex  quantity  and  is  given  by  the 
following  expression. 


*(*)=—+•'<»(*)  (5.5) 

OM 

where  a(z)  is  the  unknown  attenuation  coefficient  and  C(z)  the  sound  speed  profile. 
Since  Cb(z)=C(z)  we  write  fc(*)=Ar6(z)+ia(*).  We  also  note  that  #i(z)=/ib(z). 
Substituting  in  the  equation  (5.3),  we  obtain 


dMz) 

dz 2 


M*62M+/‘i(2)-*I2M2)=-2<V,(2M2)- 


(5.6) 


Multiply  equation  (5.6)  by  vjz)  and  equation  (5.4)  by  v(z)  and  subtract  one 
from  the  other.  We,  then,  obtain 


*>b(  *K(  z)-viz)vb"=-Mkb(  z)a{z)v(  z)vb(  z)  (5.7) 

In  the  above  equation  the  double  prime  denotes  the  second  derivative  with  respect 
to  z.  Integrating  both  sides  between  the  limits  0  to  A  leads  to  the  following 
equation. 


/  i  /"A 

vb{z)v{z)-vb{z)v(z)  \Qh  =  -2ikb{z)v(z)vb(z)a(z)dz 


(5.8) 


To  determine  u(2),*>6(2), t/(z), vb(z),  at  the  limits  2=0  and  z—h  we  make  use  of  the 
boundary  conditions,  i.e  continuity  of  pressure  and  normal  particle  velocity.  Since 
the  water  column  is  a  homogeneous  medium,  the  boundary  conditions  at  the  water 
sediment  interface  for  the  true  model  are  given  by  the  following  expressions.  R  is 
the  plane  wave  reflection  coefficient  at  the  water  sediment  interface. 

pA  0) 

fl(0)=l+/?,  F(0)  =  -1—ikJl-R)  (5.9) 

Pq 

A  similar  set  of  equations  are  obtained  for  the  boundary  condition  at  the  other 
interface.  Here  T  is  the  transmission  coefficient. 


P[h)=Teikz2k ,  P(h)=—ikz2Teikz2h 

P2 


(5.10) 


The  boundary  conditions  for  the  background  model  obtained  similarly  are  given 
below.  Rb  and  Th  are  the  reflection  and  transmission  coefficients  for  the 
background  model. 


p4(o)=i+r6,  p;(o)=-a—ikji-Rh) 

P2 


(5.11) 


The  left  hand  side  of  the  equation  (5.8)  is  now  evaluated  using  the  relationship 
between  v(z),  vb(z)  and  P(z)  and  Pb(z)  given  above  and  the  boundary  conditions. 
After  a  some  algebra  we  find  the  L.H.S.  of  equation  reduces  to  zero  at  the  upper 
limit  while  at  the  lower  limit  it  is  given  by  the  following  expression. 


(vb(z)t/(z)-v(z)vb'(z))\0=——{Rb-R) 

P  0 


(5.12: 


The  integral  equation  is  therefore, 


Po 


rh 

/  kb(z)ot{z)v{z)vb{z)dz 
Jo 


(5.13) 


This  is  a  non-linear  integral  equation  and  it  is  not  possible  to  solve  this  equation  for 
a(z)  since  v(z)  is  dependent  on  the  unknown  a(z).  As  the  attenuation  is  a  small 
perturbation  on  the  wavenumber  we  make  the  Born  approximation  v{z)=vb(z).  The 
integral  equation  then,  reduces  to  the  following. 


kJRrH)  [hHzW‘)„i 


H 

Jo 


Pb\*)dz 


(5.14) 


Pq  Jo  Pb{z) 

Since  kb(z)  and  pb(z)  are  known,  Pb(z )  is  determined  by  the  procedure  outlined  in 
Chapter  4  and  the  integral  equation  solved  for  or (z). 


The  solution  of  the  integral  equation  will  be  dealt  with  in  a  later  chapter. 
Here  we  will  continue  to  establish  the  criterion  for  the  applicability  of  the  Born 
approximation. 

5.1.2  Region  of  validity  of  the  Born  approximation 

In  studying  the  applicability  of  the  Born  approximation  we  start  with  the 
Schrodinger  type  equation  (5.3).  Putting  k{z)=kb(z)+ia{z)  we  obtain, 


d2v(z) 

dz2 


Hk2{z)+itb{z)-kx2)v(z)=-2ikba(z)v{z ) 


(5.15) 


Let  G(z,Zq)  be  the  outgoing  Green’s  function  satisfying  the  equation, 


d2G 


dz 2 


+(*,  (2l+''W_*ri  0  ^  *>2o  ^  * 


(5.18) 


Then  we  can  write, 


fh 

v(k t,z)=vb(k  ,z)+  /  2ikb(zQ)a{zQ)v(kx,z0)G(kx,z,z0)dz0 

Jo 


(5.17) 


The  total  solution  can  therefore  be  viewed  as  the  sum  of  the  solution  to  the 


homogeneous  equation  vb(z)  and  a  scattering  part  vjz). 

v(z)  =  vb(z)+ve(z) 

Using  the  above,  we  obtain 

fh 

»,(*,.*)«  /  2i*k(20)o {z0Mk  ,zQ)G(kx,z}z0)dzQ 
Jo 

fh 

/  *kb(zo)a(zo)vb(kz'Zo)G(kx'Z'Zo)dzO^ 
Jo 

fh 

/  2'kb(2o)a(zo)v9(kx’ZoWkx'Z’2o)dzO 
Jo 


(5.18) 


(5.19) 


(5.20) 


In  making  the  Born  approximation  we  write  v(*)=t;k(2).  This  implies  the  following. 


fh 

/  m(20)vt(kX’2o)G(kz'Z’Zo)dz0  <vaikx’z) 
Jo 


where, 


m(z)=2ik(b(z)a(z) 


(5.22) 


Define  the  norm  of  f(z)  as  below. 


Let  /represent  the  integral  on  the  left  hand  side  of  the  equation  (5.21).  Then, 


={/[/  rn(^a{kx,zQ)G(kx)z,zQ)dz^dz}xl2 

Jo  Jo 


But,  by  Schwartz  inequality, 


I  f  ^(zQ)v3(kx,z0)G(kx,2,20)dz0\  <  {  f  m(z0)v3(kx,z0)dzQ} 

Jo  JO 

fh 

/  G{k  ,z}zQ)dz0 
Jo 

I  f  m(2o)va(kx'zoWkz’z’zo)dzo\  ^  IMI  IKH  f  G(kx'Z'Zo)dzO 
Jo  Jo 

The  norm  of  /  therefore  satisfies  the  inequality  given  below. 


PI  <  IMI  ll»JII  f\  [  Gikt,z,z0)dz0))2dz]'l* 

Jo  Jo 

Let  G(z,zq)  have  a  maximum  value  equal  to  [G(z,z0)|ma;r •  Then, 

IMII  <  IMI  \\vt\\  h  \G(kx,z,z0)\max 

Using  this  in  equation  (5.27)  we  obtain, 

IMI  Ikll  h  \G(kx,z,zQ)\max«  ||t,,|| 
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IMI  h  \G(kx^z0)\max<1 


(5.30) 


The  condition  to  be  satisfied  for  Born  approximation  to  be  valid  is  then  given  by 


ll*ill  HI  A  |C(V.*o)L..  «' 


(5.31) 


The  applicability  of  the  Bom  approximation  depends  on  the  magnitude  of  the 
perturbation,  the  depth  of  the  sediment,  and  the  magnitude  of  the  Green’s 
function.  We  will  use  this  result  in  a  later  chapter  to  determine  the  most  suitable 
angular  aperture  for  the  inversion  scheme. 

5.2  Inverse  problem  formulation  using  the  Rytov  approximation 


We  will  in  this  section  derive  the  integral  equation  for  obtaining  the  unknown 
attenuation  coefficient  profiles  using  the  Rytov  approximation.  We  again  start  with 
the  one-dimensional  wave  equation.  For  the  sake  of  simplicity  we  have  assumed 
that  the  density  in  the  sediment  layer  is  known  and  is  a  constant. 


d2r\‘) 


+(*zW-*I2)fW=o 


(5.32) 


Let  P[z)=exp(ip(z)).  Substituting  this  in  the  above  equation  and  using  the  relation 
k{z)=kb(z)+ia{z),  we  obtain 


d?\b  dtb 

-^+{—)2-HI=b2{z)-k2)=  -  2  .'*,(*)<»(*) 
dz *  dz 


(5.33) 


Let  P^(z)  be  the  solution  to  the  wave  equation  for  the  background  model. 


-^+(W-V)A=o 

dz * 


(5.34) 


We  now  write  Pb{z)  —  exp(i/)b(z)).  Substituting  this  in  the  wave  equation  yields  the 
following. 


d\ 


+{— V+lVM-t/No 


dz£  dz 


(5.35) 


Let  ip(z)=il)b(z)+il)g(z).  Substitute  this  in  equation  (5.33)  and  then  subtract 
equation  (5.35)  from  the  result.  We  then  have  the  following  equation. 


dz  dz  dz  dz 


(5.36) 


The  left  hand  side  of  the  above  equation  is,  then,  put  into  a  more  convenient  form 
for  manipulation  as  shown  below. 


Using  equation  (5.34)  and  the  relation  Pb{z)=Pb(z)(diJ)b{z)/dz)  the  equation  (5.37)  is 


written  as  below. 


d^{ P.ib  )  d?\b  dijj.dil). 

Now  substitute  equation  (5.36)  into  the  above  to  yield 


(5.38) 


d“Pki>  dip  (z)  . 

~k,' >  =  nWK-r5— > 

dz*  dz 


(5.39) 


Multiply  equation  (5.39)  by  Pb(z)  and  equation  (5.34)  by  Pb(z)ipg(z)  and  subtract  one 


from  the  other.  We  then  integrate  over  the  depth  h  and  obtain  the  following 
equation. 


d(Phip\  dP.  L  dtp 

\"h.  pb^a  ,  Jlo  —  Pb  J  lo 


vdz 


dz 


dz 

rh  dip 


rn  dip 

=  -  /  l(— )2+2t kb(z)ot{z)]Pb2{z)dz  (5.40) 

Jo  dz 


The  left  hand  side  of  the  equation  is  evaluated  by  relating  it  to  the  boundary 
conditions  given  in  equations  (5.9), (5. 10)  and  (5.11).  The  relation  between  the  Ps 
and  ips  given  below  are  also  used  for  this  purpose. 


P(2) 

ip{z)=  In  P[z),  ipb  =  In  Pb(z),  ip  =  In  — —  (5.41) 

w 

diPg(z)  1  dP[z)  1  dPb 

-1—  - - -  (5.42) 

dz  P[z)  dz  Pft(z)  ^ z 

Using  the  above  we  find  that  the  left  hand  side  of  the  equation  vanishes  at  the 
upper  limit,  z=h.  At  the  lower  limit, 2=0,  it  is  given  by  the  following  expression. 


„  dip  ( o)  Rh~R  Ph 

F.2{z)—?—=2<kM — -  -5(X+*J2 

dz  "“(l+BKl+Rj)  ^ 


The  integral  equation  is  then 


(5.43) 


Ph  1 +Rh  M  dip 

tikdT  — =  /  K—!)2+2.'*tWaMlP62W</;  (5.44) 

p0(l+P)  Jo  dz 

So  far  in  our  derivation  no  approximations  have  been  made.  We  now  make  the 
approximation  that  ( dipjdz )2  is  small  compared  to  2kb(z)a(z).  Neglecting  this  term 
the  integral  equation  reduces  to, 


(5.45) 


k~  1 +Rb  fh  1 

—  - ~(RrR)=  /  -  kb(z)a(z)]Pb2(z)dz 

Pq  l+R  Pb 

This  equation  is  similar  to  the  equation  (5.14)  obtained  with  the  Born 
approximation  except  for  a  multiplicative  factor  in  the  left  hand  side. 

5.3  Summary 

In  this  chapter  we  have  obtained  the  linear  integral  equation  relating  the 
unknown  attenuation  profile  to  the  experimentally  determined  plane  wave 
reflection  coefficient.  We  have  applied  both  the  Born  and  the  Rytov  approximation 
to  obtain  two  similar  representation  of  the  integral  equation.  Before  we  take  up 
the  solution  of  this  integral  equation,  we  will  consider  in  the  next  chapter  the 
general  issues  involved  in  solving  this  class  of  integral  equations  and  include  in  it  a 
review  of  the  methods  in  the  literature  for  obtaining  solutions. 


have  in  effect  increased  the  magnitude  of  the  eigenvalues.  This  results  in  a  stable 
and  smooth  solution. 

To  see  how  the  addition  of  a  matrix  like  H  affects  the  eigenvalue,  we 
investigate  the  case  when  H=I  the  identity  matrix.  Consider  the  eigenvalue 
problem  forGTG.  Then, 

GTGv  =  Xv  (6.39) 

where  X  is  the  eigenvalue.  Adding  (el)v  to  both  sides,  we  obtain 

(GrG+eI)v  =  (X+eI)v  (6.40) 

The  eigenvalue  problem  for  (G  G+el)  therefore  gives  the  same  eigenvectors  but 
the  eigenvalue  has  been  increased  by  el.  The  solution  obtained  by  the 
regularisation  method  when  H=I  is  given  below. 

m  =  (G7G4-«I)“1G2cl  (6.41) 

By  eigenvector  analysis  it  can  be  shown  that  the  inverse  operator  shown  above  has 
the  following  form. 

(GrG+eI)-1  =  V„{-=-— }U„T  (6.42) 

V+e 

We  therefore  see  that  when  the  eigenvalues  are  small,  i.e  Xn<3Ce  the  diagonal  terms 
of  the  matrix  goes  as  (Xn/e),  thereby  eliminating  the  effect  of  small  eigenvalues.  On 
the  otherhand  when  Xn3>e  then  the  effect  of  e  is  ignored. 


Though  the  data  does  not  in  any  way  suggest  that  the  solution  m(x)  must  be 
smooth  function,  we  apriori  make  this  assumption.  There  could  be  other 
constraints.  But  one  selects  the  constraint  from  the  physical  nature  of  the  problem 
and  ones  guess  of  what  the  solution  should  be. 

Let  S(m(i :))  be  a  measure  of  smoothness  of  the  solution  i.e  the  smaller  the 
value  of  5,  the  smoother  m(x j  is.  Restating  the  problem  we  now  try  to  find  m(x) 
that  minimises  S(m(x))  subject  to  the  data  constraints  being  met.  A  least  squares 
procedure  can  now  be  adopted  which  minimises  |d— Gmp+\S(m(i)). 

One  measure  of  smoothness  that  is  proposed  by  TwomeyjO]  is  the  quadratic 
measure  given  below. 

rb 

S(m(i :))  =  I  (m"(x))2dx  (6.35) 

Ja 

We  can  generalise  this  by  considering  a  constraint  operator  D  acting  on  m.  We, 
then,  seek  to  minimise  the  quadratic  measure  defined  below. 

mTDTDm  =  m%m  (6.36) 

where  H=D7D.  To  obtain  a  solution  for  m  we  now  minimise  the  expression; 

mrHm+X|d-Gm|2  (6.37) 

Differentiating  with  respect  to  each  element  of  m,we  obtain  the  solution; 

(GrG+XH)m  =  GTd  (6.38) 

We  note  that  the  acceptable  solution  has  now  been  defined  as  one  which 

T 

satisfies  the  smoothness  constraint.  By  the  addition  of  the  matrix  H  to  G  G  we 


E[AmgAmgT\  =  2V / 


(6.34) 


We  note  that  the  error  in  solution  is  propotional  to  1/An2  and  the  variance  becomes 
large  when  the  eigenvalues  are  small.  Wiggins[5]  and  Jackson[6]  suggest  that  the 
effect  of  low  eigenvalues  can  be  eliminated  if  we  set  eigenvalues  lower  than  a 
threshold  level  equal  to  zero.  The  eigenvectors  associated  with  these  eigenvalues  are 
therefore  eliminated  and  the  solution  obtained  will  therefore  have  poorer  resolution. 
Eliminating  eigenvectors  corresponding  to  small  eigenvalues  is  equivalent  to 
obtaining  a  smooth  solution. 

6.2.2  Regularisation  method 

Phillips [7j,  Twomey[8j,  and  others  have  proposed  regularisation  methods  for 
obtaining  an  unique  stable  solution.  In  the  generalised  inverse  method  we  obtained 
an  unique  stable  solution  by  looking  for  a  solution  which  has  a  specific  property 
namely  that  of  minimum  norm.  The  instability  was  overcome  by  effectively 
smoothing  the  solution.  The  regularisation  method  proposed  by  Twomey[8j  is 
based  on  the  method  first  suggested  by  Phillips  [7]  and  looks  for  a  smooth  solution. 
The  measured  data  d(y)  is  available  at  discrete  points  and  therfore  it  is  defined 
only  at  these  points  and  to  within  measurement  error  as  shown  in  figure  6-1.  We 
can  therefore  say  that  d(y)  is  arbitrary  except  that  it  passes  within  the  error  bars 
associated  with  each  of  the  measurement  points.  Therefore,  there  can  be  an 
infinite  set  of  m(x)  associated  with  the  infinite  set  of  d(y).  The  ambiguity  can  only 
be  overcome  by  imposing  additional  constraints  on  m(x)  which  enables  one  to  pick 
out  of  the  large  possible  set  of  m(x)  one  that  satisfies  the  conditions  imposed  on  the 
solution.  One  such  constraint  that  can  be  imposed  is  the  smoothness  criterion. 


T 


T 


The  solution  mg  obtained  by  the  generalised  inverse  operator,  therefore,  uses  only 
the  eigenvectors  corresponding  to  the  non-zero  eigenvalues.  The  true  solution 
however  will  be  a  linear  combination  of  all  the  eigenvectors.  We  can  then  represent 
the  true  solution  by  the  following  equation. 


11* 

m  =  E  °iv.+  E  fy, 

»-l  ^n+l 

The  norm  of  the  total  solution  is  then, 


(6.30) 


m 

|m|2  =  |m(|2+  23  fi? 


(6.31) 


j=rt+l 

From  the  above  we  see  that  the  generalised  inverse  operator  yields  a  solution  with 
minimum  norm. 

Though  the  problem  of  non-uniqueness  has  been  resolved  instability  still 
exists.  Consider  the  situation  when  data  contains  error  .  Then  the  error  in 
solution  is  given  by  the  following  expression. 


Am„  =  G  "^d 

9  9 


(0.32) 


The  covariance  of  the  error  in  the  solution  can  then  be  obtained.  Here  Ej.) 
represents  the  averaging  operator. 


=  G#-1£1zldz>dT]GJ-r 


(6.33) 


Assuming  that  the  elements  of  Ad  are  statistically  independent  and  have  the  same 
variance  we  obtain, 


In  addition  there  exists  m-n  eigenvectors  Vj  such  that  Gv~=0,i=n+1, . ,m. 

Lanczos[4]  has  shown  that  the  matrix  G  can  be  decomposed  into  product  of  three 
matrices.  Using  the  decomposition  theorem  we  obtain 


G  =  IWV 


(6.26) 


Here  UN  is  an  nxn  matrix  containing  the  eigenvectors  of  GTG,  VN  is  an  mxn 
matrix  containing  the  eigenvectors  of  GTG  for  the  non-zero  eigenvalues.  A  is  a 
diagonal  matrix  containing  the  non-zero  eigenvalues.  Since  Gv—O  for  i=n+l,... 
,m,  any  linear  combination  of  these  can  be  added  to  the  solution  and  the  data  will 
still  be  satisfied.  This  is  the  cause  of  non-uniqueness  in  an  under-determined 


system. 


The  operator  defined  below  is  called  the  generalised  inverse  operator. 


G"1,  =  V^-'U  / 


Using  this  operator  .  e  obtain 


(6.27) 


=  'V1NU/<i 


(6.28) 


6.2  Solution  in  parameter  space  -  Error  free  data 


By  using  a  suitable  quadrature  scheme  the  integral  (6.3)  is  represented  by  the 
sum  given  below. 


di  =  £  "ipiTi 


(6.19) 


where  w-j  is  the  weighting  associated  with  the  quadrature  scheme.  For  the  entire 
set  of  observations  we  obtain  the  following  matrix  equation. 


d  =  Gm 


(6.20) 


where  d  is  a  nxl  vector  containing  the  observations  and  G  is  an  nxm  matrix  which 
operates  on  the  vector  m  of  dimension  mxl  representing  the  solution. 


6.2.1  Generalised  Inverse 


If  G  is  a  nxn  non-singular  matrix,  then,  the  solution  to  the  equation  is  easily 
obtained. 


m  =  G  d 


(6.21) 


When  the  integral  is  written  as  a  sum,  the  interval  Ax  is  made  small  so  that  the 
true  structure  of  the  solution  is  obtained.  This  results  in  an  under-determined 
system.  Though  the  analysis  which  follows  is  applicable  for  any  general  nxm  matrix 
we  will  consider  such  a  case  i.e.  where  n  is  less  than  m.  It  can  be  shown[4]  that  the 
square  matrices  GGT  and  GrG  define  a  coupled  eigenvector  eigenvalue  problem 


such  that 


ggV  =  x  V 


(6.22) 
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kernel  has  no  zero  eigenvalue.  One  example  of  this  is  the  Fourier  integral  which 
yields  an  unique  solution.  However,  in  practical  situations,  data  is  available  only  at 
a  limited  number  of  points.  This  situation  generally  yields  non-unique  solutions.  We 
have  to,  therfore,  deal  with  the  twin  problem  of  instability  and  non-uniqueness  in 
solution.  Many  methods  have  been  proposed  in  the  literature  for  solving  equations 
of  this  type.  All  of  them  incorporate  a  priori  information  about  the  solution  to 
overcome  the  problem  of  non-uniqueness.  The  solution  is  made  stable  by  schemes 
which  in  effect  make  the  eigenvalues  large  or  reduce  the  effect  of  the  small 
eigenvalues.  We  will  review  these  methods,  briefly,  in  the  following  paragraphs.  We 
will  also  study  how  solutions  are  obtained  when  the  data  contains  noise. 

When  data  is  available  only  at  discrete  points  the  integral  equation  takes  the 
form  in  equation  (6.3).  In  solving  the  integral  equation  two  approaches  are  possible. 
The  integral  can  be  represented  as  a  sum  by  adopting  a  suitable  quadrature  scheme 
and  then  the  solution  of  the  integral  equation  reduces  to  solving  a  set  of  linear 
equations.  A  reasonable  quadrature  scheme  can  be  adopted  if  we  have  some  prior 
knowledge  of  the  solution  from  the  physics  of  the  problem.  For  example,  if  we  know 
that  the  solution  is  smoothly  varying  we  can  assume  without  much  error  that  the 
function  is  constant  over  small  intervals  and  represent  the  solution  by  a  set  of 
parameters.  On  the  other  hand,  if  we  have  no  knowledge  about  the  solution,  such 
discretisation  can  be  erroneous  since  it  may  hide  some  structure  that  exists  in  the 
solution  unless  the  discretisation  is  made  infinitely  small.  An  alternative  is,  then, 
to  express  the  solution  as  linear  sum  of  a  set  of  basis  functions  and  determine  the 
coefficients  of  the  linear  combination.  The  class  of  methods  which  follows  the  first 
approach  we  will  call  as  providing  solution  in  parameter  space.  The  latter  we  will 
call  as  solution  in  function  space. 


»»(*)  =  2l  “JnW 


(6.15) 


Substituting  this  in  equation  (6.14)  and  using  the  orthgonality  property  of  the 
eigenfunctions  we  obtain 


%)  =  £  KanK*(y) 


(6.16) 


The  coefficients  on  are  obtained  using  the  orthgonality  property  of  the 


eigenfunctions. 


d  =  X  a  :a  =  — 

n  nn’n  . 

An 


(6.17) 


where, 


dn  =  f  d(y)<f>n{y)dy  (6-18) 

J a 

Let  d(y)  contain  an  error  e(y).  Then, 

”*(*)  =  £  “7 — K 

n=0  n 

If  the  kernel  has  small  eigenvalues  then  small  error  in  data  will  result  in  large 
errors  in  solution.  The  kernel  of  the  integral  equation  can  be  viewed  as  behaving 
like  a  smoothing  operator  which  smooths  out  wide  fluctuations  in  solution.  Such 
smoothing  operators  are  characterised  by  vhe  presence  of  low  eigenvalues.  Indeed  if 
the  kernel  is  a  delta  function  there  will  be  no  smoothing  and  a  stable  solution  is 
obtained. 

When  data  is  available  over  an  infinite  range,  unique  solution  is  possible  if  the 


*  v  '  v'  > 
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d2(y)  =  diiy)  +  N  (  G(x,y)sinux  dx 
J  O 


(6.10) 


& 


The  norm  IJ^d^d^)  is  then, 


lJidvd2)  ~  {f  l/  G(x,y)sin<jx  dx]2  dy}1/2 
J  a  J  a 

The  norm  ^(nij.mg)  is  given  by  the  foilwing  expression. 
/m(mi,m2)  =  max.  N|sinwx|  =  N 


(6.11) 


(6.12) 


For  large  values  of  u,  the  integral  (6.11)  approaches  zero  by  virtue  of  Reimann- 
Lebesque  lemma[2]  as  long  as  f*G(x,y)dx  exists.  Therefore  we  can  make  l^dvd2) 
arbitrarily  small  while  keeping  lm{ml>m2 )  large. 

To  see  how  the  kernel  of  the  integral  equation  influences  this  property,  we 
will  consider  a  special  type  of  integral  equation  for  which  we  can  obtain  solution. 
Let  us  assume  that  the  kernel  of  the  integral  equation  is  Hermitian.  Then,  by 
Hilbert-Schmit  thcorem[3]  the  kernel  can  be  represented  in  terms  of  its 
eigenfunctions.  The  eigenfuctions  are  orthogonal  and  the  eigenvalues  are  real. 


N 


G{x,y)  = 


(6.13) 


n=0 


Substituting  in  equation,  we  obtain 


fb  N 

d(y)  =  /  m(x){  (»)> dx 

J  a 


n=0 


Let, 


(6.14) 
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G{x,y)<t>{x)dx=Q 


(6.6) 


If  a  Don-trivial  solution  <f>(x)  exists  for  the  above  equation,  then  any  multiple  of  <j>(x) 
can  be  added  to  the  solution  of  the  equation  and  it  will  still  be  solution  of  the 
integral  equation.  There  can  thus  be  an  infinite  number  of  solutions  for  the  integral 
equation. 

6.1.2  Instability 

By  instability,  we  mean,  that  small  changes  in  data  can  produce  large  changes 
in  the  solution.  To  demonstrate  this  property  of  the  integral  equation  we  use  the 
proof  due  to  Tikhonov[l]. 

Consider  the  integral  equation  (6.4).  We  shall  measure  changes  in  the 
function  d(y)  by  the  norm  defined  below. 

=  I  f\dl(y)-d2(y))2dy]'l1  (6.7) 

Ja 

The  changes  in  m(x)  is  measured  by  the  following  norm 

IJm  vm2)  =  max.lmj— m2|  (6.8) 

Let  us  assume  that  the  function  mjx)  has  been  changed  to  mg(x)  given  by  the 
following  equation. 

m2(z)  =  m^zJ+TVst'^cJz)  (6.9) 


Substituting  in  the  integral  equation,  we  obtain  the  change  in  d(y). 


Linear  Inverse  Theory 


In  Chapter  5,  a  non-linear  integral  equation  was  derived  relating  the 
attenuation  coefficient  profile  to  the  reflection  coefficients.  The  equation  was 
linearised  using  the  Born  approximation  and  the  linear  equation  obtained  is  given 
below. 


rhkAz) 

W<*z>-K(*x>  =  /  -7 -a«n2(M‘'2 

Jo  oAz  1 


Jo  pb(z) 

A  similar  equation  was  obtained  using  the  Rytov  approximation.  The  unknown  in 
the  equation  is  a(z).  The  above  equation  can  be  written  as 


d(kx)  =  f  a(z)G(kx,z)dz 
Jo 


where, 


d(kT)  =  kJR^kJ-RikJ} 


Pb(z) 

This  is  a  Fredholm  integral  equation  of  the  first  kind.  In  this  chapter  we  consider 
the  issues  of  non-uniqueness  and  instability  associated  with  the  solution  of  this  type 
of  integral  equations  and  the  role  the  kernel  plays  in  this.  The  equation  given 
above  assumes  that  d(k  )  representing  data  are  available  continuously  along  a  line. 
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6.3  Solution  in  parameter  space-Data  with  errors 
6.3.1  Generalised  Inverse 

So  far  we  considered  error  free  data.  However,  in  any  measurement  there  will 
be  error.  We  can  weight  each  observation  depending  on  the  variance  of  the 
measurement.  If  o?  is  the  variance  of  the  ith  measurement,  we  weight  the  ith 
observation  by  l/<Tj.  Then  each  observation  is  related  to  the  solution  by  the 
following  equation. 


di  1 A 

-  =  “E  9iTi 


(6.43) 


We  see  that  the  data  with  the  least  variance  will  have  the  highest  weightage.  Let  S 
be  the  covariance  matrix  for  the  data.  If  the  errors  at  each  observation  are 
statistically  independent,  the  covariance  matrix  is  a  diagonal  matrix  with  the 
diagonal  terms  equal  to  1/<Tj2.  Then,  for  a  set  of  observation  we  obtain  the 
following  matrix  relationship. 


S-1/2d  =  S-1/2Gm 


(6.44) 


We  now  define  a  new  matrix  G£  =  S  1/2G  and  a  vector  dg  =  S  1/2d.  We,  then, 
obtain  the  solution  using  the  generalised  inverse  method.  If  we  have  some 
information  about  the  variance  of  the  elements  of  the  solution  vector  these  can  be 
incorporated  in  the  solution  by  multiplying  each  column  of  the  matrix  G  by  a 


weight  l/<rmj..  The  equation  to  be  solved  then  takes  the  following  form. 


d'=GW 


(6.45) 


where, 


d'  =  S_1/2d 


-121- 


G'  =  S_1/2GT1/2 
m'  =  T-1/2m 

T  is  a  diagonal  matrix  containing  <r  2  along  its  diagonal.  Using  these  relations  a 
generalised  inverse  solution  is  obtained  as  before.  Jackson [6]  has  shown  that  even  if 
the  covariance  matrices  are  not  diagonal  a  solution  can  be  obtained  which 
incorporates  the  information  on  the  covariance  of  the  data  and  the  solution. 

6.3.2  Regular isation  Method 

As  in  the  previous  method,  we  weight  each  equation  and  obtain  a  transformed 
matrix  equation.  Using  this  transformed  equation  and  the  procedure  outlined  earlier 
for  obtaining  a  regularised  solution  we  obtain[9],  the  following  solution. 

m=(GIS-1G+XH)-1GIS_1d  (6.46) 


6.3.3  Maximum  likelihood  estimatr 


The  data  measurements  are  considered  as  a  set  of  random  variables.  If  we 
make  tha  assumption  that  the  random  variables  are  distributed  in  a  Gaussian 
manner  and  S  is  the  covariance  matrix  for  these  variables,  then  the  probability 
density  function  for  the  data  is  given  by  the  following  equation. 

|S-,j1/2 


erp{— l/2(d— Gm)rS  *(d— Gm)} 


(6.47) 


For  maximising  the  probability  p(d),  we  need  to  minimise  (d-Gm)7S_,(d-Gm) 

By  a  similar  argument,  we  can  say  that  the  minimisation  of  mTm  associated  with 

0 

the  generalised  inverse  can  be  replaced  by  mT1  m.  The  solution  can  then  be 
obtained  by  the  generalised  inverse  in  transformed  co-ordinates  as  discussed  in 
earlier  paragraphs. 

6.3.4  Minimum  varaince  estimate 

In  this  approach  also  the  data  aud  the  solution  vectors  are  treated  as  a 
random  vectors.  In  the  minimun  variance  estimate  we  seek  a  linear  operator  which 
gives  an  unbiased  minimum  variance  estimate  for  the  solution.  By  Gauss- 
MarkovJlO]  theorem  the  minimum  variance  unbiased  estimate  is  given  by  the 
following  equation. 


m  =  Cm  Pid~'A 


(6.48) 


where, 


=  £l“dT| 

Cdd  =  ildd7] 

Consider  now  the  matrix  equation  given  below. 

d  =  Gm+e  (6.40) 

Here  e  represents  noise. 

/ 

Using  equation  (G.4&J  and  assuming  that  tbe  noise  and  the  solution  vector  arc 
statistically  independent  we  obtain[10] 


m  =  (GC  -1+C 

'  ee  mm 


(6.50) 


-1)GrCee-1d 

where, 

C„  =  £(eeI'l 

cmm  = 

The  stochastic  inverse  method  proposed  by  Franklin[ll]  also  yields  the  same  result. 

We  note  that  approximate  solutions  are  obtained  in  these  methods  by 
incorporating  apriori  information  about  the  statistics  of  the  solution.  Assumption 
on  the  covariance  of  the  solution  is  equivalent  to  assumptions  on  its  smoothness 
made  in  regularisation  method.  Further  addition  of  a  positive  definite  matrix  like 
C  ee*1  overcomes  the  problem  of  instability. 

6.4  Solution  in  function  space-  Error  free  data 
6.4.1  Spectral  expansion  method 

Parker(12],  arguing  that  the  inverse  problem  differs  from  parameter 
estimation  problem  in  that  the  unknown  is  a  function  with  an  infinite  number  of 
parameters,  has  provided  a  solution  in  the  function  space.  This  method  is  based  on 
the  fact  that  since  rf.  is  the  projection  of  m(x)  on  G^x)  we  can  construct  a  set  of 
orthogonal  functions  to  represent  G{ x)  and  express  m(x)  as  linear  combination  of 
these  orthogonal  functions.  To  construct  the  orthogonal  functions,  a  covariance 
matrix  f  is  formed. 


r  is  a  symmetric  matrix  and  if  the  equations  are  linearly  independent  the 
eigenvalues  are  real  and  positive,  r  is  now  expressed  as  a  product  of  three  matrices 
using  the  decomposition  theorem[4]. 

r  =  VAVT  (6.52) 

Here  the  matrix  V  contains  the  eigenvectors  and  A  is  a  diagonal  matrix  containing 
the  eigenvalues.  The  orthogonal  functions  are  then  formed  by  using  the  following 
equation. 

«■,(')  =  xr'/2Eu./'7i)  (»-«) 

j 

Expressing  the  unknown  function  as  a  linear  combination  of  the  orthogonal 
function  we  write; 

N 

m(x)  =  at4t\z)+<f>*(z)  (6.54) 

i-l 

where  $*(x)  are  orthogonal  to  0(x).  Multiplying  both  sides  by  ^(x)  and  integrating 
over  the  domain  a  to  b  we  obtain 


or, 


«,  =  \-|/2£  Vj 


(6.55) 


Having  obtained  a’s,  m(x)  is  then  evaluated  using  equation 


N 

m(x)=^2  a,4i 

1= 1 

Using  equation  (6.54)  we  now  compute  the  norm  of  m(x). 


\ 


IM*)I|2  =  f  [J2  at4t{x)+<t>*(xj\2dx 
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»=1 

or, 

||m(x)|(2  =  aar+||/||2 


(6.56) 


(6.57) 


(6.58) 


By  constructing  the  solution  using  only  <j> j  and  neglecting  the  functions  in  the 
space  orthogonal  to  it  we  obtain  the  minimum  norm  solution  as  in  the  case  of 
generalised  inverse.  As  the  a^s  are  propotional  to  Xj'1/2,  small  eigenvalues  will 
magnify  noise.  As  done  before,  we  set  the  eigenvalue  as  zero  if  its  value  is  less  than 
a  threshold  level  and  their  eigenfunctions  deleted  from  the  solution.  A  method  for 
determining  how  many  eigenfuction  are  to  be  included  in  the  solution  is  given  by 
Parker[12]. 

6.4.2  Method  of  Backus  and  Gilbert 

The  importance  that  the  kernel  of  the  integral  equation  in  obtaining  an 
unique  stable  soluion  has  been  discussed  earlier.  In  the  method  proposed  by  Backus 
and  Gilbert [13, 14],  the  kernel  of  the  equation  is  made  close  to  a  delta  function. 
Consider  the  equation  given  below. 


d~  f  GJ(x)m(x)dx  i=l,  •  ■  •  ,N 
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(6.59) 


If  by  suitable  linear  combination  of  the  kernel  we  make  it  close  to  a  delta  function 
at  some  point  xQ  then  the  value  of  the  solution  at  this  value  of  x  is  easily 
determined.  Let  a-  represent  coefficient  of  the  linear  combination.  Summing  over 


all  ap  we  obtain 
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(6.60) 
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If  A(x,xq)  is  a  delta  function  6{x— 10)  then  m(xQ)  is  equal  to  the  left  hand  side 
of  the  equation  above.  However.it  is  not  possible  to  construct  such  delta  function. 
We,  therefore,  choose  the  coefficients  ap  in  such  a  way  that  the  function  A(x,xQ) 
approaches  a  delta  function  in  a  least  squared  sense.  Various  measures  of 
’deltaness’  (spread  function)  have  been  proposed  by  the  authors[14].  For  the 
purpose  of  our  discussion  we  will  define  the  spread  function  as  given  below. 


S  =  f  {A(x,xQ)-6(x-x0)}2dx 
J  a 


(6.62) 


We  make  the  additional  assumption  that  the  function  A{x,xQ)  is  unimodular  i.e. 
J*  A(x,xQdx=l.  By  using  a  Lagrange  multiplier,  we  then  determine  the 
coefficients  which  minimises  the  spread  subject  to  the  unimodular  constraint. 
Having  determined  the  coefficients,  these  are  substituted  in  equation  ai^i 

and  the  value  of  the  solution  at  xn  determined. 


6.5  Solution  in  function  space-  Data  with  errors 

The  spectral  method  discussed  above  can  be  extended  to  take  into  account 
errors  in  data  by  weighting  each  equation.  Another  approach  which  incorporates 
the  stochastic  information  of  the  data  and  the  solution  has  been  given  by  Tarantola 
and  Nercessian  [15]  and  Tarantola  and  Valette[16|.  Tarantola  and  Nercessian[15] 
show  that  their  general  result  reduces  to  the  Backus  and  Gilbert  result  if  there  is 
no  error  in  data  and  if  we  make  the  assumption  that  no  information  exists  about 
solution. 


6.6  Non-linear  problem 

So  far  we  have  dealt  with  the  linear  problem.  The  approach  in  solving  the 
nonlinear  problem  is  to  linearise  it  by  expanding  the  kernel  around  a  base  value 
which  is  assumed  known.  We  start  with  the  non-linear  equation  given  below. 

dt  =  G[m(x),x)  (6.63) 

Now  we  expand  G^x^m^x))  around  f«0(z)  which  is  close  to  w(z)  such  that 
m(i)=m0(i)+im(i).  If  we  can  write, 


di  =  f  G{(™0(x)tx)dx+ 

Ja 

f  Di{mQ(x),x)6m(x)dx  +  0\\6m{i 
J  a 


(6.61) 


then  G^z.mjz))  is  said  to  be  Frechet  differentiable  at  m0(x).  Omitting  the  higher 


order  terms  we  can  write 


*=  *=  J  Dt{x)6m(x)dx 


(6.65) 


and  since  m^x)  is  known  ike  equation  is  now  in  linear  form  and  we  can  apply  the 
techniques  described  in  this  chapter. 

The  solution  of  non-linear  equation  starts  with  a  guess  of  the  solution  and  the 
kernel  is  expanded  about  this  and  a  correction  is  obtained.  Using  this  we  obtain  a 
new  estimate  of  the  solution  which  can  now  be  used  as  the  background.  We  note 
that  such  linearisation  is  permissible  only  if  the  kernel  is  Frechet  differentiable. 

6.7.  Summary 

In  this  chapter  we  have  reviewed  the  various  approaches  available  in  literature 
for  solving  Fredholm  integral  equation  of  the  first  kind  when  the  data  are 
available  only  at  discrete  points  and  b  contaminated  by  error.  In  the  next 
chapter  we  take  up  the  solution  of  the  integral  equation  (6.1)  and  obtain  the 
attenuation  coefficient  profile  for  the  ocean  bottom. 
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Chapter  7 

Inversion  for  Acoustic  Parameters 

In  this  chapter  we  demonstrate  the  inversion  procedure  for  obtaining  the 
attenuation  profile  for  the  ocean  bottom.  We  show  that  by  suitably  modifying  the 
scheme  corrections  to  the  sound  speed  profile  and  density  profile  for  the  bottom 
can  be  obtained  in  addition  to  the  attenuation  profile. 

First  we  establish  that  the  kernel  of  the  non-linear  integral  equation  is 
Frechet  differentiable  and  that  the  kernel  obtained  by  the  Born  approximation  is  a 
Frechet  kernel.  The  iteration  procedure  for  solving  the  non-linear  problem  is  then 
obtained. 

We  select  the  most  suitable  angular  aperture  for  input  information  based  on 
the  region  of  applicability  of  the  Born  approximation.  We  study  the  effect  that 
frequency  has  on  the  ability  to  reconstruct  so  that  this  information  can  be  used  in 
the  design  of  the  experiment.  Similarly  we  study  the  effect  of  the  magnitude  and 
the  extent  of  the  perturbation  on  the  performance  of  the  inversion  scheme.  We 
follow  this  by  demonstrating  the  reconstruction  of  different  types  of  profiles  using 
synthetically  generated  data.  A  comparison  between  results  obtained  with  Born 
and  Rytov  approximation  is  then  made. 

Reconstruction  of  the  attenuation  coefficient  profile  is  done  by  using  the 
regularisation  method  mentioned  in  Chapter  6.  We  have  not  attempted  to  carry¬ 
out  the  inversion  by  all  the  methods  described  in  Chapter  6  and  choose  one  that 
gives  the  best  results.  On  the  other  hand,  we  have  been  guided  by  the  fact  that  the 
function  we  are  trying  to  reconstruct  is  a  smooth  function.  The  regularisation 
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method  which  uses  the  smoothness  of  the  model  as  one  of  the  constraints  has  a 
strong  appeal  in  this  context  and  we  have  therefore  chosen  this  approach.  The 
results  obtained  using  this  method  are  promising. 

The  means  of  measuring  the  performance  of  the  inversion  is  developed  in  the 
context  of  the  resolving  power  theory  of  Backus  and  Gilbert[l],  Using  this  we  study 
the  resolution  of  the  inversion  performed. 

In  studying  the  question  of  errors,  we  deal  with  two  aspects  of  it.  First  we 
consider  the  situation  where  the  sound  speed  profile  and  the  density  profile  are  not 
known  exactly.  Instead  of  treating  this  difference  between  the  exact  profile  and  the 
assumed  one  as  an  error,  we  treat  it  as  an  unknown  and  include  it  as  one  of  the 
functions  to  be  determined.  We  demonstrate  this  by  simultaneously  inverting  for 
the  attenuation  profile  and  corrections  to  sound  speed  profile.  We  also  study  the 
effect  of  adding  noise  to  the  data.  Examples  of  reconstruction  using  noisy  data  are 
presented. 

7.1  Linearisation  cmd  iteration  procedure 

In  Chapter  5  we  obtained  a  non-linear  integral  equation  and  used  the  Born 
approximation  to  linearise  it.  We  will  be  using  this  linear  integral  equation  for 
performing  the  inversion.  We  now  show  that  within  the  region  of  validity  of  the 
Born  approximation  this  approximation  is  consistant  with  the  requirement  of 
Frechet  differentiability  of  the  non-linear  kernel.  To  establish  that  the  kernel  of 
the  non-linear  integral  equation  is  Frechet  differentiable  we  follow  the  analysis  of 
Parker[2].  The  integral  equation  relating  the  unknown  function  to  the  reflection 
coefficient  is  given  by  the  following  equation. 
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In  obtaining  the  above  equation  we  have  assumed  that  the  value  of  the  attenuation 
coefficient  for  the  background  model  is  zero.  We  now  write 
P[kx,z)=Pb(kx,z)+Pa{kx,z)  where  Pb{kx,z)  is  the  solution  when  a(z)=0.  Substituting 
in  the  above  equation  we  obtain, 


+  S(kI 


where 


l(z)  =  —  <*(z) 

Pb(2) 


S(kx)  =  f  'iz)Pb(kx,z)Pa(kx,z)dz 
Jo 


In  making  the  Born  approximation  we  assume  that  S(kJ  is  small  in  comparison  to 
the  first  term  and  neglect  it.  Comparing  equation  (7.2)  with  equation  (6.64)  in 
Chapter  6,  we  note  that  for  Frechet  differentiability,  S(kx)  must  be  of  the  order  of 


Define  the  norm  for  t(z)  as  given  below. 


Applying  the  Schwartz  inequality  to  S(k  )  we  find  that, 
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(7.38)  to  be  satisfied  for  larger  values  of  the  pertubation  and  sediment  thickness. 
Since  the  background  model  parameters  are  known  numerical  evaluation  of  the 
Green’s  function  can  be  performed.  From  these  computations  we  choose  the 
angular  aperture  where  the  Green’s  function  has  low  values. 

A  qualitative  feeling  for  the  acceptable  range  of  values  of  horizontal 
wavenumber  can  be  obtained  by  studying  the  simple  model  in  Figure  7-1.  For  this 
model  there  are  three  distinct  angular  region,  or  equivalently,  three  ranges  of 
horizontal  wavenumbers  that  we  will  consider.  These  are  0  <  Ar  <  k0} 
k2<  kx  <  k .,  and  fcj  <  kx  <  kQ.  For  the  first  range,  there  are  no  angles  at  which 
total  reflection  at  the  interfaces  occurs.  Energy  is  transmitted  into  the  water 
column  and  the  subbottom.  In  the  second  region  total  reflection  occurs  at  the 
bottom  interface  but  energy  is  transmitted  into  the  water  column.  In  the  region 
where  the  horizontal  wavenumber  is  greater  than  kv  there  is  total  reflection  at  the 
water  sediment  interface.  However,  the  magnitude  of  the  reflection  coefficient  at 
the  sediment/subbottom  interface  is  no  longer  unity.  Therefore  complete  trapping 
is  not  possible  for  this  particular  model.  However,  there  is  a  region  of  horizontal 
wavenumber  where  the  reflection  coefficient  magnitude  for  both  the  interfaces  is 
close  to  unity.  Near  trapping  of  energy  can  take  place  in  this  region  and  the 
Green’s  fuction  can  assume  large  values.  Figure  7-2  is  a  plot  of  the  magnitude  of 
the  Rayleigh  reflection  coefficient  at  the  two  interfaces  for  this  model.  The  region 
where  both  these  reflection  coefficients  come  close  to  unity  approximately  lies  in 
the  range  of  horizontal  wavenumbers  from  0.76  to  0.785. 

For  the  model  in  Figure  7-1,  the  Green’s  function  is  given  by  the  following 
expression. 
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Figure  7-1:  Simplified  ocean  bottom  model  to  study 
behaviour  of  Green’s  function 
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7.2.2  Angular  aperture  for  data  points 

In  Chapter  5  we  obtained  the  condition  for  the  applicability  of  the  Born 
approximation.  The  condition  to  be  met  is  given  below. 


\\kb\m\GMrnaXh<  1 


(7.38) 


The  Green’s  function  in  the  above  equation  is  for  the  assumed  background 
model.  For  any  arbitrary  sound  speed  and  density  profile  the  Green’s  function  can 
be  obtained  using  the  uniform  asymptotic  solution  for  the  wave  equation.  This 
solution  is,  however,  valid  only  when  these  parameters  vary  slowly  with  depth  and 
increase  monotonically  with  depth,  a  situation  that  is  met  in  marine  sediments.  The 
solution  for  the  Green’s  function  is  given  in  Appendix  A. 

The  Green’s  function  is  a  function  of  horizontal  wavenumber  as  well  as  the 
spatial  parameter  z  and  f,  points  representing  source  and  receiver  positions.  If  the 
Green’s  function  is  expressed  as  function  of  the  horizontal  wavenumber  there  are 
regions  in  the  horizontal  wavenumber  domain  where  the  Green’s  function  assumes 
large  values.  In  marine  sediment  the  sound  speed  at  the  sediment  water  interface  is, 
in  some  instances,  less  than  its  value  in  the  water  column.  The  sound  speed  then 
increases  monotonically  with  depth  with  a  possible  discontinuity  at  the  sediment- 
subbottom  interface.  With  this  structure,  the  sediment  layer  behaves  like  a  wave 
guide  and  complete  trapping  of  energy  is  possible.  There  will,  therefore,  be  values 
of  horizontal  wavenumber  for  which  the  Green’s  function  takes  very  large  values. 

If  we  are  to  remain  within  the  region  of  applicability  of  the  Born 
approximation  we  will  be  best  served  by  remaining  in  the  region  of  angles  where 
the  Green’s  function  has  low  values.  This  will  permit  the  inequality  in  equation 
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The  smoothness  measure  is  obtained  from  equation  (7.22)  reproduced  below. 


(7.32) 


S{a{2))  =  £K+i-2a,.+«,.-i)2 


If  we  can  have  a  measure  for  this  we  write 


(7.33) 


arHtt  <  t2 


(7.34) 


where  e 2  is  the  measure  of  smoothness.  Combining  equations  (7.32)  and  (7.33)  we 


obtain 


\Ga- d|2+— aTUa  <  2E2 

£2 


(7.35) 


We  now  find  a  solution  for  a  that  minimises  the  left  hand  side.  The  solution  is 
obtained  in  a  manner  similar  to  that  given  earlier. 


(GrG+XH)a  =  Gd 


(7.36) 


where, 


X  =  E2/c2 


(7.37) 


The  factor  E2/«2  is  a  measure  of  how  smooth  the  solution  will  be.  When  c2  is  small 
E2/c2  will  be  large  and  give  a  smooth  solution.  By  making  a  priori  assumptions 
about  E2  and  t2  we  obtain  the  Lagrange  multiplier  X. 
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decreases  in  some  region  indicating  that  the  assumption  the  the  solution  is  smooth 
over  the  entire  region  is  not  correct.  The  fact  that  such  sharp  changes  are  seen  in 
the  reconstructed  profile  even  after  making  the  smoothness  assumption  is  a  strong 
indication  that  these  features  are  not  an  artifact  of  the  inversion  scheme.  We  will, 
therefore,  be  justified  in  making  the  assumption  that  sharp  changes  do  occur  at 
these  depths  and  hence  look  for  smooth  solution  in  the  intervening  depths.  An 
example  of  this  will  be  presented  later  on. 

The  Lagrange  multiplier  X  determines  the  smoothness  of  the  solution.  The 
larger  the  value  of  X  the  smoother  the  solution.  This  means  that  larger  values  of  X 
will  lead  to  poorer  resolution.  For  lower  values  of  X  more  structure  will  be  observed 
but  at  the  same  time  the  variance  of  the  solution  will  increase.  One  approach  in 
determining  X  is  to  initially  take  large  value  and  reduce  it  with  each  iteration  as 
suggested  by  Marquadt[6].  At  each  iteration  the  residual  is  computed  and  the 
iteration  procedure  terminated  when  any  one  of  the  following  conditions  are  met. 

|d— Go|2  <  E?  (7.29) 

K+1-a„|2  <  <2  (7.30) 

Another  method  of  determining  the  Lagrange  multiplier  is  based  on  some 
prior  knowledge  of  the  upper  bound  for  the  noise  power  in  the  data  and  some 
measure  of  the  smoothness  desired  in  the  solution[7].  We  write  the  equation  to  be 
solved  in  the  following  form 

Ga  =  d+i  (7.31) 

where  s  is  the  vector  containing  the  error.  If  E2  is  the  upper  bound  for  noise 
power,  then 


to  |d— Gm|2  being  less  than  E2  (say),  the  errors  associated  with  the  numerical 
estimation  of  data.  Using  X  as  the  Lagrange  multiplier  we  seek  to  minimise 
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arHa+l/X{(d— Gaftd— Ga)} 


(7.25) 


Differentiating  with  respect  to  each  element  of  a  and  equating  it  to  zero,  we  obtain 


(G2G+XH)o  =  Grd 
a  =  (GrG+XHplGTd 


(7.26) 


(7.27) 


In  some  instances  it  may  be  necessary  to  add  additional  constraints^]  such  as 
a  minimum  norm  criterion  in  which  case  the  equation  to  be  solved  will  be  as  shown 


below. 


(G^G+XjH+XjjIJa  =  Grd 


(7.28) 


where  Xj  and  X2  are  two  Lagrange  multipliers  and  I  is  an  Identity  matrix.  In  the 
examples  considered  in  this  dissertation,  situations  needing  such  additional 
constraints  did  not  occur. 

By  using  the  method  described,  a  smooth  solution  is  obtained.  It  may  be 
noted  that  the  inversion  using  generalised  inverse  eliminates  the  instability  by 
discarding  the  small  eigenvalues.  The  small  eigenvalues  are  the  ones  that  contribute 
to  high  oscillations  in  the  solution  and  eliminating  the  small  eigenvalues  is 
equivalent  to  a  low  pass  filtering  operation  on  the  solution  to  get  a  smoothed 
version.  The  regularisation  method  also  overcomes  the  instability  by  assuming  that 
the  solution  is  in  some  sense  smooth  and  in  this  context  both  the  methods  are 
analogous. 

Having  obtained  a  solution  it  may  be  seen  that  there  are  sharp  increases/ 


Nyquist  criterion  as  the  integral  equation  looks  like  a  Fourier  transform.  Since  the 
vertical  wavenumber  varies  with  depth  and  the  angle  of  incidence,  we  will  take  the 
maximum  value  of  the  wavenumber  as  the  water  wavenumber.  Based  on  this  the 
sampling  interval  must  be  less  than  Xfl/4  where  XQ  corresponds  to  the  wavenumber 
in  the  water  column. 

The  integral  equation  is  now  written  as  a  discrete  sum 

di  =  £  (7.20) 

j 

where  w-  is  the  weighting  associated  with  the  quadrature  scheme.  For  a  set  of 
observations  the  equation  is  put  in  matrix  form, 


d  =  G  a  (7.21) 

where  the  weights  have  been  absorbed  in  G.  The  smoothness  measure  in 
equation  (7.16)  is  also  put  in  discrete  form  as  shown  below. 


5(a(*))  =  £  (a,+1~2a  •+ai._1)2 
» 

This  is  a  quadratic  measure  and  S(o(z))  is  therefore  written  as 
5(a)  =  arHa 
The  matrix  H  is  given  below. 
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The  problem  is  now  restated  as  follows.  Find  a  a  that  minimises  5(a)  subject 


synthetic  data,  this  measure  has  been  found  to  be  adequate. 


We  start  with  the  integral  equation  derived  in  Chapter  5  for  the 
reconstruction  of  the  attenuation  profile.  Since  the  reflection  coefficient  is  known 
only  for  discrete  angular  values,  the  above  equation  is  written  as 

fhkAz) 

(7.171 
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where  n  corresponds  to  the  discrete  angular  values  at  which  the  reflection  cefficents 
are  determined.  The  integral  equation  (7.17)  is  now  written  in  the  following  form 

dn=  f  a(z)Gn(z)dz  (71g) 

where 

<*„  = 

pb(z ) 

The  integral  on  the  right  hand  side  is  then  written  as  a  sum  using  a 
quadrature  scheme.  We  have  used  the  simple  Simpsons  rule  for  this  purpose.  With 
the  sediment  modelled  as  a  stack  of  homogeneous  layers  the  term  P^(kx,z)  in  the 
integrand  of  the  integral  equation  has  the  following  form  in  each  layer. 

Pt2{tx,z)  =  {Mkx)  exp(i(*i2-/:I2)1/2j)+ 

B{kx)  expt-iltj2-*/)'/2^)2  (7.19) 

Therefore,  for  the  degradation  resulting  from  converting  the  integral  into  a 
discrete  sum  to  be  small,  we  must  sample  at  points  arrived  at  by  applying  the 


small  parameter  i  as  indicated  below. 


o(r)  =  cQ0(c)+rOi(r)+  •  •  •  +  (7,15) 

We  now  successively  determine  Dq.Qj,  and  so  on.  The  convergence  of  this  scheme 
was  found  to  be  slow  compared  to  the  earlier  method. 

7.2  Reconstruction  of  the  attenuation  profile 

We,  now,  discuss  in  detail  the  reconstruction  of  the  attenuation  profile 
assuming  that  the  true  value  of  the  sound  speed  and  the  density  in  the  sediment 
layer  are  known. 

To  test  the  inversion  method  we  use  synthetic  data.  This  is  obtained  by 
solving  the  forward  problem  and  obtaining  the  plane  wave  reflection  coefficient  for 
a  set  of  angles  given  the  ocean  bottom  model.  The  data  that  will  be  used  for 
inversion  are,  therefore,  error  free. 

7.2.1  The  regularization  scheme 

Out  of  the  various  methods  described  in  Chapter  6  we  use  the  regularisation 
method.  We  assume  that  the  solution  is  in  some  sense  smooth.  Experimental  results 
available  on  the  variation  of  the  attenuation  with  depth  support  this  assumption. 
We.  therefore,  look  for  a  solution  with  the  least  structure.  The  smoothness  criterion 
is  defined  as  below [4] 

i 

fk  d:Q(:)  „ 

5(o(z))  =  /  {-^—ydz  (7.16) 
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kernel  obtained  by  Born  approximation  is  a  Frechet  kernel. 

Linearising  of  the  integral  equation  is  based  on  the  term  S(kx)  being 
negligible.  If  this  is  not  so,  then  the  solution  when  used  in  the  non-linear  integral 
equation  will  not  satisfy  data.  We,  therefore,  need  to  use  an  iterative  process.  In 
the  iterative  scheme  the  solution  obtained  in  the  previous  iteration  is  used  to 
generate  the  new  kernel.  The  integral  equation  then  takes  the  form  shown  below. 

(7.12) 


(7.13) 
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/  Hkx,<*n(*),z\dz  (7.14) 
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The  iteration  method  described  is  similar  to  the  procedure  followed  for 
solving  nonlinear  equations  by  the  Newton  method.  Convergence  of  this  scheme  is 
guaranteed  if  the  initial  guess  is  close  to  the  actual  value.  Also  the  convergence  is 
quadratic.  These  results  have  been  proved  in  the  literature[3|. 

Another  approach  is  an  iteration  scheme  based  on  higher  order  perturbation 
theory.  The  unknown  function  a(z)  is  expressed  as  a  power  series  in  terms  of  the 


fhkAz) 
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an+1U)  =  (*(z)+6ajz) 

A  variant  of  this  is  the  fixed  point  iteration.  In  this  case  we  first  write, 
&»(*)  =  a(2)-on(*) 

Substituting  this  in  equation  (7.12)  we  obtain, 


To  obtain  an  etimate  of  ||PJ|  we  asssume  that  within  the  region  of  validity  of  Born 
approximation  Pa(kx,z)  can  be  written  as  follows. 

Pa(kx,z)  =  jh^)Pb{kx,zf)G(kx,z,z!)d^  (7.6) 

J  0 

where  G(z,z’)  satisfies  the  equation; 
d^G 

—+^iHz)-kt2)G  =  (7.7) 

dz 

Applying  Schwartz  inequality  we  obtain, 


lie  II  <  IhllliniK  f (GMfdz)'!1  (7.8) 

Jo 

Using  this  we  find  that, 

l|S||  <  MI2||P»I|2{  f  (7.9) 

JO 

The  field  Pb(kx,z)  in  the  sediment  layer  is  finite.  We  will  show  in  a  later  section  of 
this  chapter  that  in  the  range  of  angles  chosen  for  input  data,  the  Green’s  function 
also  remains  finite.  The  equation  (7.9)  can,  then,  be  written  as  given  below. 


IISII 
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We  hr.ve,  therfore,  established  that  the  kernel  is  Frechet  differentiable  and  that  the 
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The  magnitude  of  the  Green’s  function  in  equation  (7.39)  is  plotted  in  Figure 
7-3  as  a  function  of  the  horizontal  wave  number.  The  source  and  receiver  positions 
are  30  m  and  10  m  respectively  from  the  water/sediment  interface.  The  acoustic 
parameters  used  in  the  computation  are  indicated  in  Figure  7-1. It  is  seen  that  the 
magnitude  of  the  Green’s  function  has  large  values  at  positions  where  the 
magnitude  of  the  plane  wave  reflection  coefficients  at  the  two  interfaces  are  close 
to  unity  as  anticipated  earlier.  The  figure  shows  also  that  in  the  region  where 
there  is  least  trapping,  i.e.  0  <kx<  k2,  the  Green’s  function  has  the  least  value. 
This  corresponds  to  pre-critical  angles.  In  the  second  region,  the  behaviour  of  the 
Green’s  function  becomes  complex  and  there  are  values  of  horizontal  wavenumber 
where  the  magnitude  of  the  Green’s  fuction  becomes  large.  When  the  waves  in  the 
sediment  becomes  inhomogeneous  the  magnitude  of  the  Green’s  function  decreases. 
The  advantage  of  using  the  information  from  the  pre-critical  range  of  angles  in  the 
inverse  problem  is  therefore  obvious.  Rays  incident  at  pre-critical  angles  will  probe 
deep  into  the  sediment  whilst  rays  incident  at  angles  above  the  critical  angle  will 
turn  within  the  sediment  layer  there  will  be  turning  points  in  the  sediment  and 
therefore  these  rays  will  have  no  information  about  medium  beyond  the  turning 
depth.  If  the  entire  range  of  angles  from  grazing  to  normal  incidence  is  used  there 
may  be  certain  angular  regions  where  the  Born  approximation  is  not  valid  and  this 
can  lead  to  degradation  in  the  reconstructed  profiles.  However,  since  an  iterative 
scheme  is  employed  the  magnitude  of  the  perturbation  gets  reduced  at  each  stage 
of  iteration  and  the  scheme  may  converge  unless  the  non-linearity  is  too  strong. 

Only  rays  in  the  precritical  region  probe  the  entire  depth  and  the  Green’s 
function  has  the  lowest  magnitudes  in  this  range.  We  have,  therefore,  used  only 
this  range  of  angles  in  the  reconstruction  of  the  attenuation  profiles. 
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Figure  7 


Magnitude  of  Green’s  function  as  function  of  horizontal 
wavenumber,  Frequency=50.0  Hz 

Source  depth=30  m,  Receiver  depth=10m 
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Figure  7-5:  Magnitude  of  Green’s  function  as  a  function  of  horizontal 
wavenumber,  Frequency—  100.0Hz 
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Source  depth=30m, Receiver  depth=10m 
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7.3  Effect  of  experimental  and  geoacoustic  parameters 
7.3.1  Frequency 

Figures  7-3  to  7-5  show  the  magnitude  of  the  Green’s  function  vs  horizontal 
wavenumber  at  different  frequencies.  In  the  pre-critical  region,  a  rough  1  /u 
dependence  is  seen.  This  result  is  not  surprising  because  in  this  region  the  effect  of 
the  interfaces  in  trapping  energy  in  least  and  the  Green’s  function  behaves  very 
much  like  the  free-space  Green’s  function  .  For  homogeneous  medium,  the  free 
space  Green’s  function  behaves  as  l/kz  where  kz  is  the  vertical  wavenumber.  This 
is  equivalent  to  \/u  dependence.  In  Appendix  A  the  expression  for  Green’s  function 
in  terms  of  the  uniform  asymptotic  solution  has  been  obtained  for  an  arbitrary 
variation  of  sound  speed  and  density.  In  the  pre-critical  range  of  angles,  there  are 
no  turning  points  in  the  sediment  layer  and  the  uniform  asymptotic  solution 
reduces  to  WKB  solution  as  shown  in  the  appendix  and  the  Green’s  function  has  an 
approximate  l/u>  dependence. 

Earlier  for  the  applicability  of  the  Born  approximation  we  noted  that  the 
inequality  to  be  satisfied  is, 


U) 


CM 


|a||*|<7(«l«’)  |  <£1 


(7.42) 


We  have  shown  that  in  the  pre-critical  angles  the  magnitude  of  the  Green’s 
function  has  an  approximate  l/u  dependence.  From  this  alone,  it  appears  that  a 
higher  frequency  is  better  since  the  magnitude  of  the  Green’s  function  is  small. 
However,  the  magnitude  of  the  perturbation  itt(z)a(z)  is  directly  propotional  to 
frequency.  Therefore  its  product  with  the  Green’s  function  is  approximately 
constant.  We  conclude,  that  the  operating  frequency  will  not  have  significant  effect 
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Figure  7-6:  Reconstruction  of  constant  attenuation  of  0.02  dB/m 
in  the  layer,  Frequency=100Hz 

1  -  True  profile,  2-  Reconstructed  profile 
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on  reconstruction.  To  see  this,  we  carry  out  inversion  using  synthetic  data 
generated  at  different  frequencies.  The  result  obtained  at  the  first  step  of  the 
iteration  process  is  at  figures  7-6  and  7-7.  Though  the  frequency  was  increased  four 
fold  it  did  not  have  any  appreciable  effect  on  the  reconstructed  profile. 


7.3.2  Attenuation  coefficient 

For  two  values  of  constant  attenuation  in  the  sediment  (0.005dB/m  and 
0.02dB/m)  inversions  were  carried  out.  The  results  are  shown  in  Figures  7-7  and 
7-8.  As  anticipated  the  larger  the  perturbation  the  poorer  the  reconstruction. 

7.3.3  Depth  of  sediment 

The  magnitude  of  the  Green’s  function  was  computed  for  different  sediment 
layer  thicknesses  while  keeping  the  other  acoustic  parameters  of  the  model  fixed. 
The  plot  of  the  Green’s  function  magnitude  as  a  function  of  horizontal  wavenumber 
is  shown  in  Figures  7-9  and  7-10.  The  Green’s  function  magnitude  remains 
approximately  the  same  for  all  the  three  cases  in  the  pre-critical  range.  Its 
behaviour  beyond  this  range  is  however  complicated.  In  equation  (7.39)  the  only 
place  where  the  sediment  thickness  enters  the  expression  is  in  the  phase  term 
exp(2»fczl/»).  Therfore  for  small  values  of  Rg  and  Rb,  which  is  approximately  true 
for  precritical  angles  the  sediment  thickness  will  have  little  effect  on  the  magnitude 
of  the  Green’s  function.  This  is  true  even  in  cases  where  the  sound  speed  and 
density  in  the  sediment  layer  varies  in  an  arbitrary  manner.  We  can,  therefore, 
conclude  that  the  sediment  thickness  has  no  significant  effect  on  the  magnitude  of 
the  Green’s  function  in  the  pre-critical  range  of  angles.  On  the  other  hand  in  the 
expression  in  equation  (7.38)  we  note  that  the  region  of  applicability  of  the  Born 
approximation  is  related  to  the  depth  of  the  sediment.  We  anticipate  that  with 
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Figure  7-8:  Reconstruction  of  constant  attenuation  of  0.005  dB/m 
in  the  layer,  Frequency  =  25Hz 


1-  True  profile,  2-  Reconstructed  profile 
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Figure  7-®:  Magnitude  of  Green's  function  as  a  function  of  horizontal 
,  wavenumber.  Depth  of  sediment—  40m 
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Figure  7-10:  Magnitude  of  Green’s  function  as  a  function  of  horizontal 
Wavenumber,  Depth  of  sediment=  120m 

Frequency*  50Hz 
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increasing  depth  the  errors  made  in  making  the  Born  approximation  increases  and 
this  will  manifest  itself  in  the  reconstructed  profiles.  Figures  7-7  and  7-11  show  the 
reconstructed  profiles  for  different  slab  thicknesses.  We  have  increased  the  slab 
thickness  four  times  leaving  every  other  parameter  the  same.  We  note  that  the 
increase  in  thickness  has  resulted  in  poorer  reconstruction. 

7.3.4  Strength  of  the  discontinuity  at  the  interface 

The  strength  of  the  discountinuity  at  the  bottom  interfaces  was  increased  by 
changing  the  values  of  the  sound  speed  in  the  sub-bottom.  The  magnitude  of  the 
Green’s  function  for  the  two  cases  are  shown  in  Figure  7-12  and  7-13.  There  is  very 
little  change  in  the  magnitude  of  the  Green’s  function  in  the  pre-critical  range  and 
therefore  the  strength  of  the  discontinuty  is  not  likely  to  have  any  appreciable 
effect  on  the  reconstructed  profiles.  However,  we  note  that  as  the  sound  speed  in 
the  sub-bottom  increases  the  angular  range  of  pre-critical  angles  decreases  and 
therefore  will  affect  the  number  of  data  points  that  can  be  used  in  the  inversion 
scheme. 

7.4  Examples  of  reconstruction 

We  now  test  the  inversion  algorithm  developed  earlier  using  different  type  of 
profiles.  The  profiles  used  are  given  in  Figure  7-14.  The  profile  shapes  have  been 
chosen  so  as  to  cover  possible  attenuation  profiles  in  marine  sediments.  The 
magnitude  of  the  attenuation  coefficient  chosen  falls  within  the  acceptable  range  of 
values  for  marine  sediments  at  the  frequency  considered. 

In  each  case  the  input  data  for  the  scheme  i.e  the  plane  wave  reflection 
coefficient,  is  computed  using  the  propagator  matrix  method  described  in  Chapter 
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Figure  7-12:  Magnitude  of  Green’s  function  as  function  of  horizontal 
wavenumber,  Sound  speed  in  sub-bottom  =  2500m/sec 


separate  regions  of  depth.  Then  the  result  in  Figure  7-20  is  obtained. 

We  now  give  an  example  of  inversion  using  one  of  the  realistic  models 
developed  for  Icelandic  Basin[9].  The  operating  frequency  in  this  instance  is  220  Hz, 
the  frequency  at  which  the  experiment  was  performed.  The  model  used  and  the 
reconstructed  profile  are  given  in  Figures  7-21  and  7-22. 

7.5  Born  and  Rytov  approximations 

The  question  about  the  applicability  of  the  Rytov  approximation  for  back 
scattering  problems  is  still  unresolved.  While  it  can  be  proved  [10]  that  the  Rytov 
approximation  is  valid  over  longer  ranges  than  the  Born  approximation  in  regions 
where  only  one  field  exists,  no  such  general  result  exists  for  regions  where  there  is 
more  than  one  wave.  According  to  Keller]  10],  where  more  than  one  wave  exists,  the 
Rytov  approximation  is  to  be  applied  to  each  wave  seperately  and  not  to  the  total 
field  if  its  validity  over  longer  ranges  than  the  Born  method  is  to  hold.  On  the 
other  hand  the  Born  approximation  can  be  applied  to  the  total  field.  This  accounts 
for  the  more  common  use  of  the  Born  approximation  for  backscattering  problems. 
However,  the  investigations  on  the  relative  merits  of  the  two  methods  for  ultra¬ 
sonic  diffraction  tomography  have  shown  that  the  Rytov  method  gives  better 
results  than  the  Born  approximation  even  when  the  Rytov  is  applied  to  the  total 
field[l  1,12].  In  the  light  of  these  comments,  we  now  carry  out  the  inversion 
without  the  iterative  scheme  and  compare  the  results  obtained  by  the  two  methods. 
The  results  are  given  in  Figures  7-23  to  7-26.  We  note  that  the  Rytov  method  is 
much  better  than  the  Born  approximation  at  shallower  depths  and  only  slightly 
worse  than  the  Born  approximation  at  deeper  depths.  This  experiment,  therefore, 
seems  to  confirm  the  observation  of  other  investigators  that  the  Rytov  method  is 
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4.  The  input  data  is  obtained  at  20  discrete  points  in  the  pre-critical  range.  For  the 
sake  of  convenience  the  input  data  is  computed  at  equal  increments  of  vertical 
wavenumber.  An  operating  frequency  of  25  Hz  was  chosen.  This  was  done  only  to 
limit  the  size  of  the  matrices  since  the  sampling  interval  is  dependent  on  the  wave 
length  in  the  water  column  as  explained  in  an  earlier  section.  The  sound  speed 
structure  in  the  sediment  was  obtained  using  the  regression  formula  given  by 
Hamilton[8].  The  density  in  the  sediment  layer  has  been  assumed  to  be  constant. 
The  values  for  these  parameters  are  given  in  Table  7-1.  The  Lagrange  multiplier 
was  suitably  chosen  and  the  iteration  continued  till  residual  became  less  than  a 
predetermined  value  based  on  the  numerical  errors  that  may  creep  into  the 
computation  of  the  reflection  coefficient. 

For  the  iteration  scheme  to  converge  the  nonlinearity  should  not  be  strong. 
Stating  differently,  the  guess  model  must  be  close  to  the  real  model.  In  Chapter  5 
we  saw  that  the  attenuation  coefficient  is  orders  of  magnitude  smaller  than  the  real 
part  of  the  wavenumber  and  this  was  the  foundation  on  which  we  based  the 
perturbation  method.  At  25  Hz  the  real  part  of  the  wavenumber  is  approximately 
0.1  m*1.  Therefore  the  initial  guess  value  for  the  attenuation  is  taken  as  zero.  For 
values  of  attenuation  assumed  the  convergence  occured  in  about  six  iterations. 

The  results  of  the  inversion  are  given  in  Figures  7-15  to  7-10.  Except  in  the 
case  of  Figure  7-19  which  is  the  recontructed  profile  for  the  discontinuous  model  at 
Figure  7-14  the  reconstruction  is  good  with  the  residuals  less  than  10  .  The  poor 
reconstruction  in  the  case  of  Figure  7-19  is  due  to  the  fact  that  the  algorithm 
assumes  that  the  function  we  are  looking  for  is  smooth  whereas  in  reality  it  is 
discontinuous.  The  points  where  such  discontinuities  are  likely  to  be,  are  seen  in 
the  Figure  7-19.  Based  on  this  evidence  we  make  the  assumption  that  the 
discontinuities  do  occur  at  these  depths  and  look  for  smooth  solutions  in  the  three 
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Figure  7-20:  Reconstruction  of  profile  in  figure  7-14(e) 
1-  True  profile,  2-  Reconstructed  profile 
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$  90 

1.7 

1616.64 

z  100 

1.7 

1628.41 

1  110 

1.7 

1640.04 

“  120 

1.7 

1651.55 

130 

1.7 

1862.92 

140 

1.7 

1674.16 

150 

1.7 

1685.28 

160 

1.7 

1696.27 

170 

1.7 

1707.14 

180 

1.7 

1717.88 

190 

1.7 

1728.50 

200 

1.7 

1739.00 

Various  profiles  adopted. 

These  are  shown  m  figure  7.14 


Sub- bottom  Density  •  19  gm/cm2 

Sound  Speed  »  1739.0  m/sec 
Attenuation  •  0.001  dB/m 


Table  7-1:  Parameters  used  for  ocean  bottom  model 
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Profile  Shops 


Description 


(a) 


(b) 


(c) 


(d) 


(e) 


OjOOS 

dB/m 


Attenuation 

Coefficient 


Constant  value  of 
attenuation*  Q005 dB/m 
in  sediment 


OOldB/m 


0  Depth  of  Sediment  200  m 


0.00675 

dB/m 

Attenuation 

Coefficient 


Linearly  decreasing  from 
.00675  dB/m  in  the  top 
1  ayer  to  .002  d  B/m  in 
.00ldB/m  the  bottom  layer 


Oepth  of  Sediment  200m 


Attenuation 

Coefficient 


.002 

dB/m 


.00675 


Linearly  increasing  profile 
from  .002  dB/m  at  the 
top  surface  to  .00675 

.OOldBAn  th«  bottom 


Depth  of  Sediment  200m 


,006dB/m 


Attenuation 

Coefficient 


jOOI 
dB/m  Q 


A  sine  curve  fitted 
between  0  and  200m  with 
a  max. of  .006  dB/m 
at  1 00m  depth 


.00ldB/m 


Oepth  of  Sediment  200m 


DOS 

dB/m 

Attenuation 

Coefficient 

.002 

dB/m 


Discontinuous  profile 
with  discontinuities 
at  40m  and  160m 


_L 


40m  160m 

Depth  of  Sediment 


200  m 


Figure  7-14:  Attenuation  coefficient  profile  shapes 
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1-  True  profile,  2-  Reconstructed  profile 


DEPTH  IN  METERS 


Figure  7-27:  Reconstruction  of  profile  in  figure  7- 14(a)  -  Rytov 

approximation 

1-  True  profile,  2-  Reconstructed  profile 
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better  than  the  Born  method  even  in  situations  where  more  than  one  field  exists. 
The  analytical  foundations  for  this  result  needs  to  be  established  to  determine  the 
limitations  in  the  use  of  Rytov’s  method.  However,  since  an  iterative  process  is 
found  necessary  to  solve  the  non-linear  problem,  both  methods  needed 
approximately  the  same  number  of  iterations  to  converge  Therefore,  neither  of  the 
two  methods  can  be  considered  to  have  a  clear  advantage  over  the  other.  The 
iterated  result  using  Rytov’s  method  for  the  profile  in  Figure  7-14  (a)  is  shown  in 
Figure  7-27 


7.6  Resolution 

Resolution  is  a  measure  of  our  ability  to  resolve  fine  structure  in  the  function 
being  reconstructed  given  that  we  have  only  a  finite  set  of  data  points.  A  method 
of  determining  the  resolvability  is  given  by  Backus  and  Gilbert,[l].  Here  the 
procedure  described  by  them  is  extended  to  cover  situations  where  the  data,  the 
kernel  and  the  function  to  be  obtained  can  all  be  complex  quantities.  The  situation 
when  the  unknown  function  is  a  complex  quantity  arises  when  we  invert  to  obtain 
the  real  and  imaginary  part  of  the  wavenumber. 

We  start  by  considering  the  linear  problem.  The  data  is  related  to  the 
unknown  through  a  linear  equation  of  the  form  given  below.  All  the  quantities  are 
complex. 


(7.43) 


Multiplying  by  a  weight  an(zQ)  and  summing  over  all  n's  we  obtain, 
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Y  dnan(2 o)  =  jf*’£  an(2o)Gn(2fidz- 


(7.44) 


or, 


d(20)=  (  ('l(z)Az,Zo)dz 
Jo 


(7.45) 


Let  d(zQ),A(z,z0)  and  -7 (2)  be  all  complex.  Equating  the  real  and  imaginary 
parts,  we  obtain 


dr=  f 

Jo 


(7.46) 


and 


fh 

di=  / 

Jo 


(7.47) 


The  subscript  r  and  i  represent  the  real  and  imaginary  parts.  If  the  coefficients 
an(2o),n~  '  *  '  >Ncan  found  such  that  A^z,z0)  is  equal  to  6(c— zQ)  and  AJ(z, zQ) 

is  zero,  then  we  will  have 


fh 

<fr(20)  =  / 

Jo 


(7.48) 


and 


^2o)=  f  '1tW6(2-zo)d2 

Jo 


(7.49) 
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<7r(20)  an<^  %(*0)  W‘H  then  he  known  exactly.  Similarly  for  each  point  we  can  find  a 
set  of  coeficents  which  will  make  Ar(z,zQ )  look  like  a  delta  function  and  AJz, zQ) 
equal  to  zero.  However,  this  cannot  be  done.  We,  therfore,  seek  to  minimise, 


z~zq)~ % 


subject  to  Ar(z,zQ )  being  unimodular  i.e. 


(7.50) 


■r(z,zo)dz=l  (7.51) 

Other  measures  of  ’deltaness’  have  been  proposed  in  literature. (1) 

From  equation  (7.44)  the  expressions  for  Ar(z,zQ)  and  A((z,z0)  are  obtained. 
Substituting  for  Ar(z,zQ)  and  At{z,zQ)  in  the  equation  (7.46),  we  minimise, 

JQ  E  an^GnM)  +  JQ  E  an(2o)Gn,(2)l2d2 

n  n 

+X  f  E  (7  52) 

Jo 

n 

where  X  is  the  lagrange  multiplier.  The  derivative  of  this  expression  with  respect  to 
each  of  the  coefficient  aQ  is  set  to  zero.  A  set  of  n  linear  equations  is  obtained 
which  are  then  solved  to  obtain  the  coefficients.  After  some  algebra  the  following 
equation  is  obtained. 

2(*r+*,)a+Xu  =  2Gr(z0)  (7.5o) 


where 


'Vnm  Jq  nr'  ’  mr'  > 


Wnm  =  f  Gni^Gmx<z)dz 
JO 

a={al(*o)>  ‘  ‘  *  'an(zoftT 
GMo)  =  {Glr("o)-  •  •  *  >G„r(*0)}r 


Then, 

2tfra+Xu  =  2G  r(z0) 
Or, 


a  =  #-'(Gr(20)-W2l 


We  now  use  the  constraint  relation  given  below 


aru=l 


Sustituting  for  a,  wc  obtain 


uT#1Gr(20)-Xu/2l  =  1 


(7.58) 


UT^>GJ20)-1 

-! :■■■ }«| 

U1*  *U 


(7.59) 


Having  found  an,  we  can  construcL4r(z,z0)  and  A^z^Zq).  Ar(z,z0)  is  a  measure 
of  the  resolution.  Further  At{z, zQ)  is  a  measure  of  the  effect  of  the  imaginary  part 
on  the  recontruction  of  the  real  part  and  vice-versa  and  is  therefore  a  measure  of 
contamination. 

The  above  method  is  applicable  when  the  equation  is  linear.  In  the  non-linear 
problem,  the  resolution  kernels  are  obtained  by  the  method  suggested  by 
Parker[13j.  The  linearisation  assumption  is  based  on  the  last  term  in  equation  (7.2) 
being  negligible.  If  this  is  not  small,  the  solution  we  obtain  will  not  satisfy  data.  By 
iterating  we  obtain  an  acceptable  solution  and  when  such  a  solution  is  obtained  we 
can  say  that  the  term  0||7||2  is  negligible.  The  Frechet  kernels  obtained  with  this 
solution  is  the  used  to  generate  the  resolution  kernels. 

7.6.1  Measure  of  resolution 

Resolution  length  defined  as  given  below  is  a  measure  of  resolution(14j. 


Rl\z{ 


1  rhA2(z,zQ ) 

hJo  A2(z0,zq) 


(7.60) 


While  resolution  length  at  each  point  z  gives  a  measure  of  local  resolution  one  can 
also  obtain  a  measure  for  global  resolution. 


v-s 


1  fh 

GRL=-  I  Rl^z^dzy  (7.61) 

hJo 

The  global  resolution  measure  can  be  used  to  determine  the  effect  of  various 
parameters  in  inversion. 

For  the  various  examples  of  the  reconstruction  of  the  attenuation  coefficient 
profile  the  resolution  kernels  have  been  constructed  adopting  a  procedure  similar  to 
the  one  discussed  above  except  that  the  function  7 (z)  is  real.  There  will,  therefore, 
no  contamination  from  the  imaginary  part  i.e  we  need  to  determine  only  A^z^z^). 
These  results  are  given  in  Figures  7-28  to  7-31.  The  resolution  lengths  for  the  cases 
are  given  in  Table  7-II. 

The  average  resolution  length  is  approximately  .07  which  corresponds  to  14 
m.  The  layer  thickness  has  been  taken  as  10m  and  within  each  layer  the 
attenuation  is  taken  as  constant.  A  resolution  length  of  14  m  can  be  considered  as 
adequate  in  resolving  the  features  of  the  model. 

We  will  now  see  how  the  aperture  size  and  the  number  of  data  points  affect 
resolution.  We  will  use  the  global  resolution  measure  to  study  this  effect. 

7.6.2  Variation  of  aperture  size 

Keeping  the  number  of  data  points  same,  the  aperture  size  is  varied  and  for 
each  case  the  resolution  length  and  the  global  resolution  length  are  caiculated.  The 
resolution  length  and  global  resolution  measure  are  given  in  Table  7-111  .  As  the 
angular  aperture  increases  the  resolution  improves.  With  small  angular  aperture, 
the  separation  between  data  points  is  small.  In  the  ray  picture  the  rays  probing  (he 
medium  will  be  close  to  each  other  leading  to  near  dependency  of  the  rows  of  the 
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matrix  we  are  trying  to  invert.  This  illconditioning  has  to  be  overcome  by 
discarding  eigenvectors  corresponding  to  low  eigenvalues  and  therefore  leads  to 
poor  resolution.  However  as  the  aperture  size  increases  the  near  dependency  of  the 
rows  gets  reduced  and  the  resolution  improves. 

7.6.3  Variation  of  the  number  of  data  points 

For  this  experiment,  the  resolution  length  and  global  resolution  length  were 
obtained  for  number  of  data  equal  to  11,17  and  20.  These  are  given  in  Table  7-IV. 
As  the  number  of  data  points  increase  the  resolution  improves.  This  is  what  one 
would  anticipate  as  the  small  number  of  points  makes  the  system  underdetermined 
with  a  number  of  eigenvalues  close  to  zero.  A  smooth  solution  is  obtained  by 
utilising  a  relatively  small  number  of  eigenvectors  which  degrades  the  resolution. 

7.7  Errors 

So  far  we  have  assumed  that  the  sound  speed  and  density  for  the  sediment 
are  known  exactly.  In  the  sequence  of  processing  of  the  measured  data  as  shown  in 
Figure  (3-7)  of  Chapter  3,  the  sound  speed  and  density  values  are  obtained  using 
direct  inverse  algorithms.  These  algorithms  assume  that  the  sediment  is  lossless  and 
since  this  assumption  is  not  true,  it  will  give  rise  to  errors  in  the  reconstructed 
values  of  the  sound  speed  and  density.  To  show  this  we  do  the  following  synthetic 
experiment.  Merab[15]  has  developed  an  inverse  scheme  based  on  the  Gelfand- 
Levitan  method  for  obtaining  the  sound  speed  and  density  profiles  for  the  sediment. 
We  generated  reflection  coefficients  synthetically  using  a  lossy  model  and  this  was 
used  in  the  inversion  scheme  to  obtain  the  value  of  sound  speed.  The  density  in  the 
sediment  was  assumed  known.  The  procedure  was  repeated  with  a  different  set  of 


Exact  Profile 


input  data  corresponding  to  the  lossless  case  with  every  other  parameter  in  the 
model  remaining  the  same.  The  lossy  and  lossless  models  used  are  shown  in  Figure 
7-32.  The  results  of  inversion  is  given  in  Figures  7-33  and  7-34.  Even  with  the 
attenuation  coefficient  as  small  as  0.001  dB/m  there  is  observable  error  in  the 
reconstructed  value  of  the  sound  speed.  However,  the  result  obtained  can  be  used 
as  the  initial  guess  model  for  the  sound  speed.  We  will  call  such  errors  in  our 
knowledge  of  sound  speed  and  density  in  the  sediment  as  model  errors  and  show 
how  such  errors  can  be  dealt  with. 

The  second  form  of  errors  are  those  in  the  data.  The  data,  i.e.,  the  plane- 
wave  reflection  coefficients  are  not  measured  directly  but  are  obtained  by  suitably 
processing  some  other  measured  quantity.  For  the  scheme  proposed  in  Chapter  3 
the  directly  measured  quantity  is  the  pressure  field  and  this  is  then  processed. 
Errors  that  occur  in  the  measurement  of  this  quantity  will  therefore  be  propagated 
and  will  finally  appear  as  error  in  the  reflection  coefficient.  We  will  call  these  data 
errors  and  show  by  way  of  examples  how  these  affect  the  reconstruction  of  the 
attenuation  profile. 

7.7.1  Errors  in  the  model 

To  show  the  effect  of  errors  in  the  sound  speed  and  density  we  use  the  model 
in  Table  7-1  .  The  input  data  is  generated  using  the  correct  value  of  the  sound 
speed.  However  while  performing  the  inversion  a  slightly  different  valu?  of  sound 
speed  is  assumed.  A  similar  test  is  done  for  error  in  density  as  well.  The 
degradation  that  occurs  in  the  solution  is  shown  in  Figures  7-35  and  7-36.  Table 
7-V  gives  the  true  sound  speed  and  the  value  assumed  for  reconstruction.  In  the 
case  of  density  the  constant  value  of  1.7gm/cc  in  the  sediment  layer  is  changed  1.75 
gm/cc.  In  both  cases  the  attenuation  profile  used  for  recontruction  is  in  Figure 
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Figure  7-38:  Effect,  of  error  in  density  on  reconstruction 
1-  True  profile,  2-  Reconstructed  profile 


Consider  the  case  where  there  is  an  error  in  the  background  model  for  the 
sound  speed  profile.  We  can,  then,  relate  the  wavenumber  of  the  true  model  to  the 
wavenumber  of  the  guess  model  by  the  expression  given  below. 

AC(z) 

C^z) 

k(z)—kb+6k(z)  (7.62) 

Since  6k(z)  is  not  known  it  cannot  be  subtracted  out.  An  alternative  approach  is  to 
treat  this  as  an  unknown  and  invert  to  obtain  Sk(z)+ia(z). 

We  now  take  up  the  general  case  where  we  consider  errors  in  the  sound  speed 
and  density  profiles  of  the  guess  model.  The  parameters  for  the  true  model  and  the 
background  model  are  now  related  by  the  following  equation. 


k(z)  —  /rfr(2)+M:(2)+ta(2) 

(7.63) 

P(z)  =  Pb(z)+0(z) 

(7.64) 

Our  aim  is  to  obtain  an  integral  equation  which  on  inversion  will  yield  the 
corrections  to  the  sound  speed  and  density  profiles  in  addition  to  the  attenuation 
coefficent  profile.  To  do  this  we  make  the  assumption  that  the  density  profile  and 
its  first  derivative  are  continuous  across  the  interfaces.  We,  then,  use  the  procedure 
developed  in  Chapter  5  to  derive  the  integral  equation  for  this  case.  Using 
equations  (5.3)  and  (5.4),  we  express  fi( z)  as  a  sum  of  /^(z)  and  a  perturbation  of  it 
as  given  by  the  following  equation. 


k(z)= - [1+ 

Cb(z) 


p(z)  =  Vb(z)+6(z) 


(7.65) 
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where, 


d2/3(z)  dpb(z )  dp(z)  d2pb{z)  dp(z) 


=  —  +  3_  —  -|- 


+  3- 


dz 


-W) 


(7.66) 


The  integral  equation  then  is, 


=  fh^[^b(z)(6k(z)+ia{z))+6(z)) 

JO  pb{z) 

Pb2(kx,z)dz  (7.67) 

This  equation  can  be  solved  to  obtain  2kb(z)6k(z)+6(z)  and  2kb(z)o(z).  To  seperate 
5k(z)  and  S(z)  we  use  the  method  due  to  Coen(14).  The  experiment  is  performed  at 
two  frequencies  and  for  each  case  the  integral  equation  is  solved.  Let  the  solutions 
obtained  be  t]l  and  t}2  respectively.  Then  the  real  part  of  the  solutions  can  be 
expressed  as  follows. 

u2 

»?1  M)  =  2-11r~AC(z)+6{z)  (7.68) 

C?(z) 

O 

f/2r(2)  =  22—AC(z)+6{z)  (7.69) 

cb3M 

Subtracting  equation  (7.69)  from  (7.68),  we  obtain 


AC(z)  = 


CfeV)  1ir(z)-12r(Z) 


2  2 
CJ.  — < 

"1  w2 


(7.70) 


We  obtain  2$C(z)  and  6(z)  from  the  above.  Having  obtained  6(z)  we  then  obtain  /?( z) 
by  using  the  equation  (7.66).  The  boundary  condtions  for  the  problem  are  /?(0)=0 
and  £’(0)=0. 


Now  consider  the  case  when  there  is  only  error  in  the  sound  speed.  In 
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functions  being  reconstructed.  We,  therefore  ,  seek  to  minimise  the  sum  of  the 
smoothness  measures  subject  to  the  data  constraints  being  met.  The  solution  then 
takes  the  form 

[GTG+X(Hj+H2)]m  =  Grd  (7.75) 

where  Hj  and  H2  will  depend  oh  the  smoothness  measure  for  the  functions. 

Using  such  a  scheme,  simultaneous  reconstruction  of  the  sound  speed  and 
attenuation  profile  was  done.  Tables  7- VI  and  7- VII  give  the  details  of  the  exact 
and  initial  guess  values  for  the  two  cases  considered.  The  exact  and  reconstructed 
profiles  are  in  Figures  7-37  to  7-40. 

While  considering  the  effect  of  frequency  on  the  inversion  for  the  attenuation 
profile  we  stated  that  the  frequency  has  no  effect  on  inversion.  However,  when 
simultaneous  inversion  is  performed  the  perturbation  associated  with  it  is  the 
magnitude  of  5k(z)+ia(z).  Since  5k  is  propotional  to  the  square  of  frequency,  the 
product  5kG(z,z’)  will  be  propotional  to  frequency.  This  leads  us  to  conclude  that  in 
this  case  a  lower  frquency  is  better  for  reconstruction. 

The  case  where  there  are  errors  in  our  knowledge  of  the  subbottom 
parameters  is  treated  in  Appendix  B 

7.7.2  Errors  in  data 

We  have  already  mentioned  that  to  model  the  error  in  the  data  we  will  need 
to  know  the  proceesing  that  has  been  done  to  obtain  the  data,  namely,  the  plane 
wave  reflection  coefficients.  In  the  scheme  envisaged  in  Chapter  3  the  first  stage  is 
to  obtain  the  Green’s  function  by  carrying  out  a  Hankel  transform  operation  of  the 
pressure  field.  Mook[16]  studied  the  effect  of  adding  stationary  white  gaussian  noise 
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Figure  7-37:  Simultaneous  reconstruction  -  sound  speed  profile;  case  (i) 
3-  True  profile,  4-  Reconstructed  profile 
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Table  7- VII:  True  and  background  model  parameters  -  case(n) 


tea&aasift'  GSomotk  i&a&ssfet  aMBSflffi  «aas^  jaiaaafe 


Exact 

initial 

Exact 

Soundspeed  m/s 

Guess  value  m's 

dB/m 

1517.56 

1511.00 

0.005 

1530.39 

1523.00 

0.005 

1543.14 

1535.00 

0.005 

1555.74 

1547.00 

0.005 

1568.20 

1559.00 

0.005 

1580.52 

1571.00 

0.005 

1592.69 

1583.00 

0.005 

1604.74 

1S9S.00 

0.005 

1616.64 

1607.00 

0.005 

1628.41 

1619.00 

0.005 

1640.40 

1631.00 

0.005 

1651.55 

1643.00 

0.005 

1662.92 

1655.00 

0.005 

1674.16 

1667.00 

0.005 

1685.28 

1679.00 

0.005 

1696.27 

1691.00 

0.005 

1707.14 

1703.00 

0.005 

1717.88 

1715.00 

0.005 

1728.50 

1727.00 

0  005 

1739.00 

1739.00 

0.005 

Initial  guess 
dB/m 


Table  7- VI:  True  and  background  model  parameters  -  case(i) 


Mragfr? 


equation  (7.65)  we  set  $(2)=0.  We  then  obtain  the  integral  equation  given  below. 


fhkb[z) 

ikM{R(kx)-Rb(kx)}  «  /  -^-{6k(z)+ia(z)}Pb2(kx,z)dz 

JO  Pb{z) 


(7.71) 


For  discrete  values  of  kx  we  obtain, 


dn~  [  m(z)Gn(z)dz  (7-72) 

Jo 

where 

*n  =  =  Wn+*K)„ 

G„W=GJ*>+iGmW 

m(z)=mr(2)+im((z) 

with  Gnr(z)  and  Gni(z)  being  the  real  and  imaginary  parts  of  Pb2(kn,z).  If  the  data 
are  available  at  discrete  data  points,  we  we  employ  a  quadrature  scheme  and  obtain 
a  matrix  equation  of  the  form  given  below. 

d=Gm  (7.73) 

where  d  is  the  vector  (dr,dj]T,  m  is  the  vector  [mr,mj]T  and  G  is  the  matrix  given 
below. 


G= 


(7.74) 


To  solve  this  matrix  equation  we  again  use  the  regularisation  procedure  and 
impose  constraints  on  the  solution  based  on  smoothness  criteria  for  the  two 
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to  the  pressure  field  .  If  the  noise  is  zero  mean,  the  expected  value  of  the  output  of 
the  Hankel  transform  will  not  be  corrupted  since  the  Hankel  transform  operation  is 
a  linear  one.  However  the  variance  is  found  to  be  non-st  at  ion  ary  and  the  noise 
power  is  concentrated  around  kf  =0  where  k  is  the  horizontal  wavenumber.  The 
non  stationarity  of  the  noise  is  simply  that  it  is  dependent  on  kr  value  and  this 
makes  it  difficult  to  model  noise  in  the  reflection  coefficient.  However  we  make  the 
assumption  that  within  the  small  angular  aperture  that  we  use  for  the  data  that 
the  noise  is  stationary  and  gaussian.  With  this  assumption  we  define  signal  to 
noise  ratio  as  given  below. 


SNR=\0log.J~) 


(7.76) 


For  reflection  coefficients  at  a  set  of  N  discrete  points  we  have 


(7.77) 

(7.78) 


Using  this  measure,  errors  corresponding  to  40  dB  and  50  dB  was  added  to  each 
data  point  and  inversion  of  the  attenuation  profile  carried  for  diiferent  type  of 
profiles.  The  results  obtained  for  the  different  cases  are  shown  in  Figures  7-41  to 
7-46.  This  experiment  was  performed  to  assess  whether  the  inversion  method  is 
stable  in  the  presence  of  noise.  We  note  that  for  the  noise  levels  used,  the  scheme  is 
stable  and  yields  reasonable  results. 


In  this  chapter  we  demonstrated  that  the  linear  integral  equation  derived  in 
Chapter  5  can  be  solved  to  obtain  the  attenuation  coefficient  profile  for  the  ocean 
bottom.  Though  we  started  with  the  problem  of  requiring  to  determine  the 
imaginary  part  of  the  wavenumber,  a  perturbation  of  the  real  part,  we  showed  that 
the  scheme  can  yield  perturbations  in  both  the  real  and  imaginary  part  thereby 
enabling  determination  of  all  the  three  acoustic  parameters.  We  demonstrated  this 
by  determining  the  attenuation  coefficient  profile  and  the  sound  speed  profile.  We 
studied  the  effect  of  experimental  and  acoustic  parameters  on  inversion.  The 
resolving  power  theory  of  Backus  and  Gilbertfl]  was  used  to  assess  the  resolution 
obtainable  with  the  data  and  examined  the  effect  of  the  angular  aperture  and  the 
number  of  data  points  on  inversion.  Finally  we  showed  by  examples  that  the 
inversion  method  is  stable  in  the  presence  of  noise. 

In  the  next  chapter,  we  examine  briefly  another  approach  to  obtain  the 
acoustic  parameters  of  the  ocean  bottom  which  is  suitable  in  the  shallow  water 


context. 
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Chapter  8 

Acoustic  Parameters  from  Eigenvalues 

In  this  chapter  we  will  show  that  in  the  shallow  water  context  a  first  order 
perturbation  theory  can  be  used  to  connect  perturbations  in  the  modal  eigenvalues 
to  the  acoustical  properties  of  the  ocean  bottom.  Such  a  formulation  leads  to  an 
integral  equation  of  the  type  we  have  considered  so  far  in  earlier  chapters  and 
therefore  the  procedures  recommended  for  the  solution  of  such  equations  can  be 
used  here  as  well. 

Frisk  and  Lynch[l]  have  proposed  that  an  experiment  similar  to  the  one 
described  in  Chapter  3  can  be  carried  out  in  the  shallow  water  and  measurement  of 
the  pressure  field  taken.  The  experimental  configuration  is  shown  in  Figure  8-1  .  If 
the  pressure  field  obtained  is  now  Hankel-transformed  to  obtain  the  Green’s 
function,  distinct  peaks  can  be  observed  at  the  modal  eigenvalues.  This  has  been 
demonstrated  by  Frisk  and  Lynch[2]  using  real  data.  Figure  8-2  is  the  Greens 
function  computed  by  carrying  out  the  Hankel  transform  of  the  pressure  field 
measured  in  an  experiment  conducted  off  the  Nantucket  sound.  The  position  of  the 
poles  of  the  Green’s  function  are  easily  seen.  These  occur  at  the  modal  eigenvalues. 

If  the  medium  is  absorbing  the  modal  eigenvalues  will  become  complex  and 
the  real  and  imagimary  parts  of  the  eigenvalues  need  to  be  found.  Frisk  and 
Lynch(l]  give  a  method  by  which  this  can  be  done.  We,  therefore,  assume  that  from 
the  measurements  of  the  pressure  field  an  estimate  of  the  modal  eigenvalues  have 
been  obtained.  The  eigenvalues  of  the  modes  depend  on  the  acoustic  parameters  of 
the  bottom  and  therefore  have  this  information.  The  problem  then  is  to  develop  a 
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Figure  8-1:  Experimental  configuration  for  shallow  water  experiment 
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Figure  8-2:  Hankel  transform  of  measured  pressure  field 


procedure  to  reconstruct  the  geoacoustic  parameters  of  the  bottom  from  a 
knowledge  of  the  modal  eigenvalues. 

8.1  Derivation  of  integral  equation 

Consider  the  model  shown  in  Figure  8-3.  In  the  shallow  water  environment  we 
assume  that  the  sound  speed  in  the  water  column  is  constant  and  all  the  acoustical 
parameters  in  the  water  column  are  known.  We  are  now  required  to  determine  the 
density  .sound  speed  and  attenuation  profiles  for  the  bottom. 

We  assume  that  from  archival  data  or  other  sources  of  information  we  are 
able  to  get  approximate  values  for  these  parameters.  This  is  used  as  the  initial 
guess  model  for  the  bottom  and  a  perturbative  approach  used  to  construct  an 
integral  equation  which  will  relate  the  changes  in  the  eigenvalues  between  that  of 
the  guess  model  and  the  real  model  to  the  acoustic  parameters  of  the  bottom. 

Making  the  assumptions  of  circular  symmetry  to  a  point  source  field  and 
horizontal  stratification,  the  equations  for  the  normal  mode  is  given  by  the 
following  equation. 

d2UJz)  X  dUlz) 

—f—  +/>(*)<— )'-r 1 — +(fc2(^)-fc„2)l/„(r)=0  (8.1| 

dz *  p(z)  dz 

where  k(z)  is  the  complex  wave  number  for  the  medium,  i.e  absorption  in  the 
medium  is  accounted  for  by  making  the  wave  number  complex.  kn  is  the  modal 
eigenvalue  and  Un  the  modal  function  for  the  nth  mode.  Making  the  substitution 
vn(2)=P~l^{2Wn(2)  we  obtain, 


d\(z) 


+{k2(z)+»(z)-kn2)vn(z)= 0. 


vtViA  vCVCVCV.''-' 


>  J> 


Let  Cb(z)  and  pb(z)  be  the  initial  guess  model  for  the  ocean  bottom.  The 
attenuation  is  assumed  to  be  zero  for  the  guess  model.  Then, 


Hkbz(z)+Pb(z)-kz2)vnb(z)=0 


dzl 


where, 


Pb1,2A'W  „ 

=  T~W 


Let, 


**)  =  Vnb(Z)  +  w  J2) 


and 


(8.3) 


(8.4) 


kiz)=kb+ia(z);  n{z)=pb(z)+6{z)  (8.5) 

Substituting  the  above  into  equation  (8.2), we  obtain 

wn6W+wn/+lA:62(2:)-2fc6^)^2^d)^2)+,'2^0^r^  + 

A*6+^^)-*n2)t;n^)==0  (8-6) 


Adding  and  subtracting  knb 


vn(z )  to  the  left  hand  side  of  equation  (8.6),  we  obtain 


Vn6W+Vn/+lA:62+m6(2)-fcn62)Vn(2)+ 

(2A:6(z){2g(d)^2)+iQ(2)}+5(2)]t;n(z)-{fcn2-*62}i;n(2)=0 


(8.7) 


Subtracting  equation  (8.7)  from  equation  (8.3),  we  obtain 


vn"+\kbW+»bW-knbKM) 

=-{2kb(6k{z)+ia(z))+S{z)}vn{z)+(kn2-knb2)vn{z) 
Taking  the  complex  conjugate  of  equation  (8.3),  we  obtain 


d\h  (*) 


—f — +(kb2(z)+»b (2)~kx2Kb  (2)=°  (8  °) 

dz 

Multiplying  equation  (8.8)  by  v nb*  and  equation  (8.9)  by  vng  and  subtracting  we 
obtain 


+(kn2~knb2K(2Kb(2) 


(8.10) 


We  now  integrate  over  the  entire  depth.  The  left  hand  side  then  becomes  equal  to 
zero.  The  resulting  equation  is  then, 


( kn2~knb2)f  vn(zKb  *(z)dz= 
JQ 


f  \2kb(z){6k{z)+ia{z)}+%z)}\{z)vnbWdz 

Jo 


(8.11) 


Using  the  first  order  approximation,  vn(z)=unfc(2),  we  obtain 


*UJZ)  Unb  (z) 


WTO 


Since  Un s  are  orthogonal,  we  obtain 


k  2-k  2 

Kn  Knb 


■=  /  —  [2*,(2){«(2)+m(!:))+ii(2)||l/n,(2)|2rfZ 
Jo  p(z) 


(8.13) 


*«=*„6+**=*n6  +(»Jr+‘'Wi 


(8.14) 


Then, 


(8.15) 


Substituting  these  relations  in  equation  (8.13)  we  obtain  the  following  integral 
equation  for  the  unknown  parameters. 


kJSkJr  rl 


-tk^S^+mzWPnfdz 


(8.16) 


knh(6kn)i  f0®1  o 

”  =  /  —Uz)ak(z)\Unfdz  (8.17) 

Nn  Jo  i>bW 

If  the  experiment  is  performed  at  two  frequencies  6k{z)  and  6(<r)  can  be  seperated. 


!=  r'—kp)ck{z)\UJ2dz 
Jo  pAz) 


The  approach  is  the  same  as  the  one  shown  in  Chapter  7.  We  can  recover  the 
corrections  to  density,  sound  speed  and  attenuation  in  this  manner.  The  similarity 
between  the  equations  (8.16)  and  (8.17)  and  those  used  in  connection  with  the  deep 
water  problem  is  evident.  Therefore  the  methods  suggested  in  Chapter  6  can  be 
used  to  solve  the  above  equations  and  obtain  the  geo-acoustic  parameters  of  the 
bottom. 


8.2  Inversion  for  sound  speed  using  synthetic  data 

Unlike  the  problem  where  the  reflection  coefficent  at  many  angles  of  incidence 
are  available  as  data,  in  this  case  the  number  of  data  points  are  limited  by  the 
number  of  modes.  Further  the  modes  have  a  decaying  field  in  the  bottom  and 
therefore  the  penetration  of  the  modes  is  small.  The  changes  in  the  eigenvalues  are 
affected  by  the  sediment  properties  only  in  the  top  layers.  If  the  source  frequency  is 
reduced  then  the  modes  penetrate  deeper  but  the  number  of  modes  gets  reduced, 
thereby  reducing  the  number  of  data  points  available  for  inversion. 

To  test  the  validity  of  this  method,  the  synthetic  data  were  generated.  A 
model  for  the  ocean  bottom  was  assumed  and  the  modal  eigenvalues  determined 
numerically  by  using  a  search  routine  which  searches  for  the  location  of  the  poles. 
For  the  guess  model  the  modal  eigenvalues  were  simiiarly  obtained.  The  bottom 
was  assumed  to  be  lossless.  Having  thus  obtained  the  difference  in  the  modal 
eigenvalues  inversion  algorithm  was  used  to  reconstruct  the  bottom  acoustic 
parameter [3].  In  this  case  it  was  assumed  that  the  density  in  the  bottom  was 
known  leaving  only  the  bottom  sound  speed  to  be  determined. 

To  study  the  effect  of  the  magnitude  of  the  perturbation  two  cases  were 
studied;  one  in  which  the  perturbation  was  taken  as  50  m/s  and  the  other  in  which 
it  was  100  m/s.  The  exact  and  guess  models  are  given  in  Figures  8-4  and  8-5. 

The  modal  eigen  values  were  computed  at  three  different  frequencies  50 
Hz,  100  Hz  and  200  Hz.  The  data  from  the  two  frequencies  50  Hz  and  100  Hz  were 
then  seperately  used  in  an  inversion  algorithm.  The  results  obtained  are  given  in 
Figures  8-6  and  8-7.  In  order  to  increase  the  number  of  data  points,  the  data  from 
all  three  frequencies  were  combined  and  used  in  the  inversion  procedure  In  this 
case  we  further  made  the  assumption  that  the  location  of  the  sub-bottom  is  known. 


solutions.  When  Q(z)  is  a  smoothly  varying  function  and  if  there  is  only  one  turning 
point,  the  solution  can  be  written  in  termr  of  the  uniform  asymptotic  solutions. 
Let  g^z)  and  g0(z)  be  the  two  linearly  independent  solutions.  Then, 


9x(z)  =  Q  ll\z)uxl\z)Ai{-u) 
9 2iz)  =  Q~ll2(z)vXl\z)Bi(-v) 


(A3) 


(A4) 


where, 


^M3/2a,e(2)}2/3 


=  [zQm 

Jzt 

zt  is  the  turning  point  i.e.  Q(zt)= 0.  Ai  and  Bi  are  the  Airy  functions. 

We  now  write  the  solutions  in  the  two  regions  0  <  z  <  z?  and  z1  <  z  <  h  in 
terms  of  these  solutions. 


<7(5, 2/)  =  Agl(z)+Bg2{z)  0  <  2  <  z1 


(AS) 


G(z,z!)  =  Cgl(z)+Dg2(z)  z1  <  z  <  h 


(A.6) 


These  solutions  must  satisfy  the  boundary  conditions  at  2=0  and  z—h.  Further  the 
function  G{z,J)  must  be  continuous  at  z=z!.  dG/dz,  however,  is  discontinuous  at 


In  the  region  2  <  0,  the  region  is  homogeneous.  The  field  in  this  region  is  an 
outgoing  wave  of  the  form  TQ  exp(»^z02).  kiQ  is  the  vertical  wave  number  in  the 
region  2  <  0.  Applying  boundary  conditions  of  continuity  of  pressure  and  normal 
particle  velocity  at  2=0, 
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Appendix  A 


Consider  the  model  in  Figure  (2-16).  We  now  derive  the  Green’s  function 
which  satisfies  the  following  equation. 

d2G 

—  +  (k2(z)+p(z)-k2)G  =  Siz-z*)  0  <  z,z*  <  h  (Al) 

dz“ 

where, 
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u> 
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Px!\z)  p\z ) 
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Equation  (A.l)  can  then  be  written  as  below. 


d2G 


dz 


-  +  uj2Q2(z)G  =  qz-z*)  0  <  z,zf  <  h 


(A2) 


where 


=  { 


/i(z)  sin  2»0 


C2(z) 


w 


C  2 
co 
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To  solve  for  the  Green’s  function,  we  note  that  the  Green’s  function  satisfies  the 
homogeneous  form  of  equation  (A. 2)  when  z  Therefore  in  the  regions  where 

z  2*.  we  write  the  solution  as  the  linear  combination  of  the  two  independent 


0.1.2  Effect  of  shear 


In  the  development  of  the  reconstruction  scheme,  we  assumed  the  sediment  to 
be  a  fluid.  However,  the  sediment  supports  shear  and  part  of  the  congressional 
wave  energy  is  converted  to  shear.  This  will  manifest  itself  as  additional  loss.  The 
effect  of  shear  on  the  reconstruction  therefore  needs  investigation.  Also  possibility 
of  extending  the  method  to  give  information  about  shear  is  another  aspect  to  be 
investigated. 

0.1.3  Shallow  water  problem 

Only  very  preliminary  investigation  has  been  done  on  the  method  proposed 
for  reconstruction  of  acoustic  parameters  from  modal  eigenvalues.  Further  work 
needs  to  be  done  in  this  regard. 

0.1.4  Inversion  from  the  point  source  pressure  field 

The  input  information  for  the  inversion  scheme  is  the  plane  wave  reflection 
coefficent  while  the  measured  quatity  is  the  complex  point  source  pressure  field. 
The  measured  pressure  field  is  processed  to  get  the  input  data.  Errors  in  input  data 
as  well  as  assumptions  made  in  the  processing  of  the  data  will  affect  the  quality  of 
the  reflection  coefficents  obtained.  For  example  it  has  been  shown  that  one  the 
sources  of  error  in  the  Hankel  transform  operation  is  the  variation  in  source  height. 
A  scheme  which  takes  the  measured  pressure  fields  as  data  will  be  less  prone  to 
such  processing  errors  and  therefore  the  possibility  of  developing  such  a  scheme 
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part  of  the  solution  affects  the  reconstruction  of  the  real  part  and  vice-versa. 

We  showed  that  the  reconstruction  of  the  attenuation  coefficent  is  sensitive 
to  errors  in  the  sound  speed  and  density  and  then  developed  and  demonstrated  a 
method  where  the  errors  are  treated  as  additional  unknowns  in  the  inversion. 
However,  when  reconstruction  of  the  sound  speed  is  performed,  the  errors  increase 
with  frequency.  Hence  the  accuracy  of  the  reconstruction  will  be  limited  by  the 
frequency  at  which  the  experiment  is  performed. 

The  r' feet  of  additive  noise  in  the  input  data  was  investigated  and  in  the  case 
of  the  reconstruction  of  the  attenuation  profile  with  noisy  data,  we  noted  that  the 
scheme  is  stable. 

In  Chapter  8,  we  derived  the  integral  equation  for  obtaing  the  acoustic 
parameters  of  the  ocean  bottom  from  modal  eigenvalues.  This  method  is  applicable 
to  shallow  water  wave  guides.  We  demonstrated  the  inversion  using  synthetic  data 
and  one  set  of  data  acquired  in  a  field  experiment.  The  preliminary  results  show 
promise  and  further  work  is  continuing. 

The  foundations  laid  in  this  thesis  provide  a  basis  for  future  work  in  other 
areas  some  of  which  are  presented  in  the  next  section. 

0.1  Suggestions  for  future  work 

0.1.1  Testing  the  validity  with  real  data 

We  used  primarily  synthetic  data  in  this  thesis.  The  validity  of  the  scheme 
needs  to  be  tested  fully  using  real  data. 


problem  using  a  uniform  asymptotic  solution  was  also  indicated  in  the  chap' or. 

In  Chapter  5,  we  derived  the  linear  Fredholm  integral  equation  of  the  first 
kind  for  obtaining  the  attenuation  coefficent  profile  for  the  sediment  layer.  The 
linear  equation  was  obtained  using  the  Born  and  Rytov  approximations.  We 
derived  an  expression  for  the  region  of  validity  of  the  Born  approximation  and 
showed  that  it  is  dependent  on  the  magnitude  and  extent  of  the  perturbation.  A 
review  of  the  methods  available  for  the  solution  of  this  class  of  integral  equation 
was  made  in  Chapter  6  and  an  iterative  method  for  solving  non-linear  integral 
equation  indicated. 

In  Chapter  7  we  established  that  the  kernel  of  the  non-linear  integral 
equation  is  Frechet  differentiable  and  that  the  Kernel  obtained  by  the  Born 
approximation  is  a  Frechet  kernel.  We  then  arrived  at  the  optimum  aperture  angle 
for  input  information  by  studying  the  region  of  validity  of  the  Born  approximation. 
We  noted  that  this  corresponds  to  pre-critical  angles.  We  also  showed  that  if  the 
sound  speed  and  density  are  known  exactly,  and  only  the  attenuation  coefficient  is 
desired,  then  the  frequency  at  which  the  experiment  is  conducted  does  not  affect 
the  reconstruction. 

The  inversion  procedure  to  obtain  the  attenuation  coefficent  profile  was 
demonstrated  using  synthetic  data  for  different  types  of  profiles  including  a 
discontinuous  one.  The  solution  was  obtained  by  making  an  a  priori  assumption  on 
the  smoothness  of  the  solution.  The  ability  of  the  data  to  resolve  fine  details  of  the 
solution  was  studied  by  computing  the  Dirichlet  resolving  kernel  using  the  resolving 
power  theory  of  Backus  and  Gilbert.  In  the  case  of  noise  free  data  the  method 
yielded  adequate  resolution.  We  also  derived  the  resolving  kernel  when  the 
unknown  is  a  complex  quantity  and  obtained  an  expression  for  the  resolving  kernel 
and  a  contamination  term.  The  contamination  term  determines  how  the  imaginary 


Chapter  9 
Conclusions 


In  this  thesis,  we  developed  a  perturbation  method  for  obtaining  the  profiles 
of  the  compressional  wave  speed,  the  compressional  wave  attenuation  and  density 
in  marine  sediments.  We  inferred  these  parameters  using  an  inverse  procedure  that 
has  as  its  input  the  plane  wave  reflection  coefficent  as  function  of  angle  of  incidence 
at  a  fixed  frequency.  The  Born  approximation  was  applied  to  obtain  a  linear 
Fredholm  integral  equation  of  the  first  kind.  A  numerical  method  for  solving  the 
integral  equation  using  a  priori  information  on  the  smoothness  of  the  solution  was 
adopted  and  the  inversion  procedure  tested  with  noisy  and  noise  free  synthetic 
data.  In  both  cases  the  scheme  yielded  stable  and  acceptable  results. 

In  Chapter  2  we  studied  the  propagation  of  plane  waves  in  marine  sediments 
and  developed  an  acceptable  ocean  bottom  model.  We  noted  that  no  acceptable 
theoretical  model  as  yet  exists  to  describe  the  variation  of  the  attenuation 
coefficient  at  low  frequencies  in  marine  sediments.  Further,  we  noted  that 
dispersion  of  the  propagating  waves  can  be  substantial  for  coarse  sediments.  This 
led  us  to  propose  a  field  experiment  using  a  mono-chromatic  source  for  determining 
the  acoustic  properties  of  the  ocean  bottom.  In  Chapter  3  we  briefly  described  the 
experiment  and  enumerated  the  processing  steps  to  be  used  to  obtain  the  acoustic 
parameters. 

In  Chapter  4  we  dealt  with  the  forward  problem  and  obtained  &  numerically 
stable  propagator  matrix  algorithm  to  compute  the  plane  wave  reflection  coefficient 
and  the  wave  field  in  the  sediment.  An  alternate  method  of  solving  the  forward 


-237- 


8.4  References 


1.  G.V.  Frisk  end  J.F.  Lynch,  "Shallow  water  waveguide  characterisation 
using  Hankel  transform,"  J.  Acoust.  Soc.  Am., 76,  203-216  (1084). 

2.  G.V.  Frisk,  J.F.  Lynch,  and  J.A.  Doutt,  "Shallow  water  waveguide 
characterisation  using  Hankel  transform  -  Preliminary  results,” 
J.  Acoust.  Soc.  Am.,  76(sl)S(94),  (1684). 

3.  S.D.  Rajan,  M.S.  Wengrovitz,  J.F.  Lynch,  and  G.V.  Frisk,  "A 
linear  inverse  scheme  for  the  extraction  of  ocean  bottom 
geoacoustic  parameters  from  modal  eigenvalues,”  J.  Acoust.  Soc. 

Am.,  76(Sl)  S(64),  (1084). 

4.  E.L.  Hamilton,  "  Geoacoustic  modelling  of  the  sea  floor,  ”  J.  Acoust. 
Soc.  Am.,  68(5),  1313-1340  (1080). 


WATER  (13.9m,  ISOVELOCITY) 


Figure  8-10:  Inversion  of  real  data 
1-  Initial  guess,  2-  Reconstructs  -ofile 
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The  inversion  method  was  similar  to  the  one  described  in  the  last  chapter,  namely 
that  we  looked  for  a  smooth  solution.  The  results  are  given  in  Figures  8-8and  8-9. 

The  following  observations  are  made  on  the  basis  of  results  obtained, 


1.  In  the  case  of  high  frequency,  reconstruction  of  only  the  top  layer  is 
possible.  At  lower  frequencies  reconstruction  to  a  deeper  depth  can  be 
obtained. 

2.  When  data  obtained  at  three  different  frequencies  are  used  as  input 
data  the  situation  improved  considerably  as  more  data  points  are 
available. 

3.  If  we  add  more  information  to  the  scheme  such  as  the  speed  in  the 
subbottom  and  the  location  of  this  interface,  the  results  of  inversion 
are  better. 


8.3  Inversion  of  real  data 

The  data  obtained  in  the  field  experiment  was  processed  and  the  modal  eigen 
values  determined  at  two  frequencies  140  Hz  and  220  Hz.  Inversion  for  the  sound 
speed  structure  of  the  bottom  was  carried  out  using  the  method  proposed  earlier. 
The  result  obtained  is  shown  in  Figure  8-10.  The  profile  for  the  bottom  is  similar  to 
the  one  proposed  by  HamiIton(4]  for  sandy  sediments. 

Work  on  this  approach  to  obtain  the  acoustic  parameters  of  the  bottom  is 
continuing  and  a  more  detailed  investigation  will  be  necessary  before  any 
conclusions  can  be  made  on  the  efficacy  of  the  proposed  method. 
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Figure  8-0:  Reconstructed  profile  using  data  from  all 

three  frequencies-  case(ii) 


1-  Initial  guess,  2-  Reconstructed  profile, 
3-  Reconstructed  profile  after  one  iteration 
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Figure  8-8:  Reconstructed  profile  using  data  from  all 

three  frequencies-  case(i) 


i-  Initial  guess,  2-  Reconstructed  profile, 
3-  Reconstructed  profile  after  one  iteration 
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Figure  8-7:  Reconstructed  profile  -  case(ii) 
frequencies-  50Hz  and  100Hz 

Initial  guess,  2-  Reconstructed  profile  using 
100  Hz  data,  3-  Reconstructed  profile  using  50Hz  data 
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Figure  8-6:  Reconstructed  profile  -  case(i) 
Frequencies  -  50Hz  and  100Hz 


1-  Initial  guess,  2-  Reconstructed  profile  using 
100  Hz  data,  3-  Reconstructed  profile  using  50Hz  data 
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^(0)  +  Bg2(0)  =  TQ 


(A.  7) 


—  {Ap1/2{z)gl{z)+  Bpl/2{z)g2(z)Y\0  =  —  TQ 
P(0)  P0 


(AS) 


We  now  assume  that  p(z)  and  its  first  derivative  are  continuous  across  the 
interface.  With  this  assumption,  we  obtain 


B  p0g^)  +»*i0/»1/2(o)ff1(o) 

( AM) 

We  now  write  the  solution  in  the  region  z  <  2*  as  follows. 

<7(2,*')  =  B{Rsgx(z)  +  g2{z)} 

(AAO) 

where  Rg  =  A/ B  and  is  given  by  the  equation  (A.9). 

boundary  conditions  at  z  =  h,  we  obtain  for  2  >  2*, 

Similarly  by  applying  the 

G(‘S)  =  C{},(*)  +  R^z)} 

(All) 

where 

P09i(h)-ikz2pl,2{h)gi(h) 

’■YV' 


Sragtra 


i 


5.;: 


i>l(z)  =  Rt9\(2)  +  92(2) 


i>2(z)  =  9\(z)  +  Rb92iz) 


(A14) 


(A  15) 


Cj  is  a  constant  and  is  equal  to  1  /yV{rpl(z),tl)2(z)}.  The  Wronksian  W{tl)l(z),rp2(z)  is 


given  by 


=  V>2(  *#/(*)  - 

Substituting  for  ip^z)  and  if>2(z)  fr°m  equations  (A. 14)  and  (A.15),  we  obtain 


(A. 16) 
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W[ip  ^2)^(2))= - —— 

u  l~RsRb 

The  solution  for  the  Green’s  function  is  therefore, 
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G(z,z>)= - 
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G(z,z>)= - 
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(A- 18) 


(A- 19) 


If  there  are  no  turning  points,  then  the  uniform  asymptotic  solution  reduces  to  the 
WKB  solution.  Then, 
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z  <  z1 


(A20) 


(A.21) 


With  this  solution  for  g^z)  and  g^z),  we  can  show  that  Rg  and  Rb  look  like  the 
Rayleigh  reflection  coefficents.  Let  us  consider  Rg. 


-s " ' 


(A. 22) 


0l'(O)  = 


Q1,2(  0) 
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For  a  smoothly  varying  medium  {Q'(z)lQ\z)}<&  1  and  in  the  frequency  range  we 
are  considering  this  term  is  negligible.  Therfore, 

0l'(°)=*M(°)}1/2  (>4.23) 

Similarly, 

?2'(°)=-i{wQ(°)}1/2  (A24) 

But  wQ(0)=^(0)  where  kg(0)=k2(0)—kz2.  Using  these  relations  in  equation  (A.2) 
we  obtain, 

A  Pokz(°)-P(0)kto 

—  = -  (A. 25) 

B  P0kz(0)+p(0)kM 

Compare  this  with  the  equation  for  the  Rayleigh  reflection  coefficent  R  between 
two  medium  with  densities  pQ  and  px  and  sound  speeds  CQ  and  Cy  For  an  incident 
wave  in  the  medium  with  parameters  pl  and  Cj 
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Appendix  B 


In  this  appendix  we  indicate  how  errors  in  the  acoustic  parameters  of  the 
subbottom  can  be  dealt  with.  In  this  we  follow  the  method  of  Cocn[l]. 

Consider  the  case  where  the  error  is  in  the  value  of  the  attenuation  coefficient 
in  the  subbottom.  We  express  the  true  value  of  the  attenuation  coefficent  as  a  sum 
of  a  background  value  and  a  small  error  term. 

a2  =  a2b+Sa  (B.l) 

Proceeding  in  exactly  the  same  manner  as  in  Chapter  5,  an  integral  equation  is 
obtained. 
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For  discrete  values  of  angle  we  write 


dn  =  j  a(z)Pb n(2)dz  + 

JO 

f  8aPbn(z)dz 

J  h 


(B.2) 


(B.3) 


Since  the  unknown  6a  is  a  constant  it  is  taken  out  of  the  integral  and  the  following 
equation  obtained. 
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a 
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where 


dn=  f  a(2)Pbn^dz+SaIn 
JO 


(BA) 


K  = 


(B.5) 


At  some  other  angle  we  write 


dm=  f  «(z)Pbm2(z)dz+6aIn 
Jo 


(B.  6) 


Eliminating  6a  from  these  two  equations  the  following  relation  is  obtained. 


(B.7) 


For  each  data  point  l,....,n,  the  contribution  from  the  second  term  is  subtracted 
out  by  using  the  (n+l)th  data  point  as  shown  above.  The  equations  are  now  solved 
for  a(z)  as  before.  Having  obtained  a(z),  6a  is  then  computed.  Error  in  other 
parameters  can  be  treated  in  a  similar  manner. 
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